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Abstract 

The  Air  Force's  Prompt  Global  Reach  concept  describes  the  desire  to  have  a 
capability  to  reach  any  target  within  a  9000  nautical  mile  distance  within  two  hours  of 
launch.  To  meet  this  objective,  much  effort  is  being  devoted  to  hypersonics  and  re-entry 
vehicles.  A  hypersonic  vehicle  typically  has  poor  aerodynamic  efficiency,  and  therefore, 
can  only  make  small  alterations  to  its  trajectory.  Given  the  limited  maneuverability, 
computational  modeling  is  used  to  generate  trajectories  before  launch  to  strike  intended 
targets.  In  addition  to  endpoint  (target)  constraints,  waypoints  may  be  necessary  for 
reconnaissance,  multiple  payloads,  or  other  logistical  activities.  Once  the  optimal 
trajectory  is  solved  to  satisfy  the  endpoint  and  waypoint  constraints,  the  next  question 
asked  is,  "Where  else  can  the  vehicle  go  while  still  meeting  the  mission  objectives,  and 
what  is  the  penalty  for  making  such  maneuvers?"  The  result  of  this  research  is  a  direct 
numerical  solution  technique  for  mapping  the  sensitivity  of  the  terminal  state  as  a 
function  of  additional  waypoint  location  while  satisfying  vehicle  dynamics,  control 
limitations,  and  all  previously  defined  waypoint  constraints.  Multiple  cases  are 
presented  including  a  simple  endpoint-to-endpoint  scenario  and  a  waypoint  included 
scenario,  with  a  Gauss  pseudospectral  solver  as  the  direct  numerical  solver. 
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Optimal  Re-entry  Trajectory  Terminal  State  Due  to  Variations  in  Waypoint 

Locations 

I.  Introduction 

1.1  Motivation 

The  United  States  Air  Force's  Concept  of  Operations  (CONOPS)  includes  Global 
Strike  (GS)  and  Global  Persistent  Attack  (GPA).  These  concepts  are  the  result  of  an  Air 
Combat  Command  study  conducted  in  2002,  in  which  the  United  States  Air  Force 
(USAF)  queried  the  industry  about  which  technologies  and  research  would  be  needed  to 
have  a  Long  Range  Global  Precision  Engagement  (LRGPE).  "The  key  focus  area  is  on 
the  capability  to  strike  targets  anywhere  and  anytime  within  the  global 
battlespace. .."[1],  The  key  focus  was  further  refined  to  include  the  element  of  time. 
Therefore,  research  is  being  dedicated  to  hypersonic  and  re-entry  vehicles  capable  of 
engaging  any  target  on  the  Earth  within  12  hrs,  i.e.  Prompt  Global  Strike. 

A  more  recent  progress  report  on  where  the  United  States  stands  in  terms  of 
Prompt  Global  Strike  was  given  by  General  Chilton  in  his  statement  to  the  House 
Armed  Services  Committee  on  27  February  2008  [14],  He  states  that  we  have  a  "Prompt 
global  strike  capability  on  alert  today,  but  it  is  configured  only  with  nuclear  weapons." 
Nuclear  weapons  have  a  large  range  of  destruction  that  may  include  assets  we  wish  to 
protect,  thereby  precluding  this  option  as  a  means  of  deterrence.  General  Chilton  goes 
on  to  say,  "The  capability  we  lack  is  the  means  to  deliver  prompt,  precise,  conventional 
kinetic  effects  at  inter-continental  ranges.  The  ability  to  hold  at  risk  sites  in  otherwise 
denied  territory  is  a  key  element  of  our  strategic  deterrent  capability."  This  is  a  gap  in 
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the  United  States  deterrent  capability  that  is  receiving  high  levels  of  interest  and 
attention  for  a  near-term  solution  [14], 

The  motivation  highlighted  above  addresses  the  terminal  objective  of  GS  and 
GPA,  but  another  facet  of  the  problem  involves  optimal  trajectory  generation.  The 
targets  are  considered  to  be  time  critical,  therefore  the  wish  is  to  minimize  the  time  of 
flight  to  target.  However,  there  are  constraints  to  be  considered  such  as  political 
boundaries  and  borders,  physical  obstacles  like  mountains,  and  the  desire  to  fly  over 
specific  waypoints  for  multiple  payload  departure  points,  communication  relay,  and 
other  tasks.  Coupled  with  these  tasks  are  the  limitations  of  the  vehicle:  vehicle 
dynamics,  heating  constraints,  g  load,  etc. 

In  addition  to  the  military  needs  cited  above,  NASA's  Integrated  Space 
Transportation  Plan/Program  (ISTP)  seeks  to  expand  the  civil  and  commercial  reach 
into  space  for  the  coming  decades  [5].  This  plan  emphasizes  Reusable  Launch  Vehicles 
(RLV).  Presently,  there  are  two  general  thoughts  for  retrieving  RLVs,  either  recover  the 
RLV  downrange  of  the  launch  site  and  transport  it  back  to  the  launch  facility  or 
fly/  glide  the  RLV  back  to  the  launch  site.  The  latter  option  is  preferred  since  it  does  not 
require  off-site  resources  to  recover  the  vehicle,  and  the  vehicle  can  otherwise  land  and 
be  serviced  immediately  for  a  quick  turn-around  time. 

1.2  Problem  Description 

The  long-term  goal  of  this  research  is  to  provide  a  tool  to  mission  strategists 
capable  of  providing  near  real-time  feedback  of  all  feasible  waypoints  given  an  initial 
condition  and  final  target.  Such  a  scenario  faced  by  a  strategist  would  include 
identifying  a  finite  number  of  targets  or  waypoints,  determining  if  the  waypoints  are 
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feasible,  and  evaluating  the  associated  costs.  The  choice  of  one  waypoint  will  reduce  the 
subset  of  further  feasible  waypoints,  and  this  process  will  continue  until  the  subset  of 
remaining  feasible  waypoints  is  zero.  A  trajectory  minimizing  the  time  to  final  target 
will  determine  how  the  RV  maneuvers  between  discrete  waypoints. 

In  hypersonics,  it  is  somewhat  commonly  known  that  small  deviations  off  the 
nominal  trajectory  can  throw  the  vehicle  significantly  off  course.  The  same  is  true  if  the 
trajectory  is  purposely  modified  to  move  a  waypoint  or  avoid  a  no-fly  zone.  Therefore, 
the  type  of  problem  to  be  addressed  herein  is  a  sensitivity  analysis  to  waypoint  location. 
Although  closely  related  to  perturbation  and  uncertainty  theory,  sensitivity  theory  is 
interested  in  characterizing  changes  in  some  parameter  with  respect  to  a  specific  event. 

In  particular,  the  effects  on  terminal  state  with  respect  to  changes  in  waypoint  location 
will  be  investigated. 

In  summary,  the  goals  of  this  research  are  to: 

Use  a  high  fidelity  simulation  of  a  hypersonic  re-entry  vehicle  with  a  fixed  initial  condition  and 
preselected  final  target  position  to: 

a.  Determine  a  minimum  time  flight  trajectory  satisfying  all  imposed  vehicle  constraints. 

b.  Investigate  the  sensitivity  on  the  primary  mission  objective  (time  to  target)  as  the 

number  and  location  of  intermediary  waypoints  are  imposed  on  the  trajectory. 

In  addition  to  the  above  goals,  this  research  intends  to  increase  the  model  fidelity 
for  a  more  realistic  simulation.  In  general,  simplified  assumptions  in  the  model 
dynamics  are  justified  by  the  need  for  comparatively  short  computational  time.  The 
decision  to  make  simplifications  to  the  model  dynamics  must  be  balanced  with  accuracy 
and  minimizing  error.  To  illustrate,  imagine  a  compass  onboard  a  ship  in  the  open 
ocean  that  has  an  accuracy  of  ±T .  As  the  ship  approaches  a  course  heading  change,  the 
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magnitude  of  the  heading  change  only  needs  to  be  known  to  the  nearest  degree  since  the 
compass  is  only  good  to  ±1° .  However  if  the  ship  must  navigate  through  narrow 
channels  or  risk  running  aground,  then  both  the  compass  accuracy  and  required 
heading  change  need  to  be  known  with  greater  precision.  Two  different  types  of 
compasses  could  be  used  for  each  of  the  applications;  i.e.,  one  compass  that  uses  the 
magnetic  field  lines  with  a  correction  based  on  location  for  a  quick  but  not  as  precise 
heading,  while  the  other  compass  uses  a  GPS  system  which  requires  a  few  seconds  to 
calculate  the  heading  based  on  the  difference  between  the  ship's  present  position,  to  its 
position  ten  seconds  prior.  Hypersonic  vehicles  can  be  likened  to  the  ship  navigating 
through  narrow  channels,  and  requiring,  therefore,  precise  navigation.  The  previous 
work  upon  which  this  research  is  built  used  a  simplified  dynamics  model  that,  while 
unrealistic,  nevertheless  proved  useful  in  analysis  of  an  optimal  trajectory.  Therefore 
the  second  goal  of  this  research  is  to: 

Apply  high  fidelity  3-D  dynamics,  the  results  of  which  may  he  compared  with  results  derived 
using  the  simplified  model. 

The  challenge  to  finding  an  optimal  re-entry  trajectory  in  general  is  further 
complicated  by  vehicle  control  limitations,  material  properties,  and  human  limitations. 
Such  restrictions  are  applied  to  ensure  a  survivable  trajectory  for  both  the  vehicle  and 
crew,  if  manned.  Figure  1  below  shows  the  flight  corridor  for  the  space  shuttle  with  two 
different  models.  The  graphic  on  the  left  is  from  the  flight  test  program,  and  on  the  right 
is  the  operational  angle-of-attack  profile.  Although  the  pre-computed  trajectory  may 
satisfy  the  computer  model  dynamics,  it  may  violate  the  operational  constraints. 
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Accounting  for  all  dynamics  and  constraints  in  trajectory  generation  becomes  important 
for  vehicle  (and  crew)  survival. 


DRAG 

ACCELERATION. 

FT/SEC2 


FLIGHT  CORRIDOR 


Figure  1.  Shuttle  Entry  -  Angle-of- Attack  Profiles  [18]. 
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II.  Previous  Research 


This  research  extends  the  research  carried  out  by  Maj.  Jorris  [23].  His  research 
provided  a  "synergistic  solution  encompassing  three  key  technologies;  trajectory 
generation,  waypoint  satisfaction,  and  threat  or  no-fly  zone  avoidance."  His  solution 
used  a  Gauss  Pseudospectral  Method  (GPM)  form  of  dynamic  optimization.  His 
technique  lays  out  the  groundwork  for  solving  the  optimal  trajectory,  but  it  is  only  part 
of  the  research  required  to  tackle  this  new  problem.  The  other  key  area  of  research 
addresses  the  solution  sensitivity  to  changes  in  mission  parameters  and  incorporates  a 
spherical  Earth  model  and  its  rotation. 

2.1  Dynamic  Optimization 

Dynamic  optimization  has  the  goal  of  finding  the  optimal  control  variables  such 
that  an  objective  function,  J,  is  minimized  or  maximized  without  violating  state 
dynamics,  boundary  conditions,  and  equality  and/ or  inequality  path  constraints.  The 
general  optimal  control  problem  is  posed  below. 

Minimize 

J  =  <I)(x(f0),f0,Jc(f/),f/)  +  |j  '  L(x(t),u(t),t)dt  (2.1) 

subject  to  the  dynamic  constraints 

x(t)  =  f(x(t),u(t),t)  (2.2) 

the  event  constraints  (boundary  conditions) 

</>(x{t0),t0,x(tf),tf)  =  0  (2.3) 

and  the  path  constraints 

C[x(t),u(t),t)<  0  (2.4) 
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where 


xe  C(M") 

^er  (2.5) 

C  g  c(M"') 

and  n,  p,  and  m  are  the  number  of  states,  event  constraints  and  path  constraints 
respectively. 

This  type  of  problem  has  been  addressed  from  many  angles.  For  simple 
problems,  analytical  solutions  have  been  obtained  and  for  others  only  numerical 
solutions  exist.  Numerical  methods  for  solving  the  optimal  control  problem  are  solved 
via  indirect  or  direct  methods.  The  indirect  method  is  characterized  by  formulating  the 
first-order  optimality  conditions  from  the  objective  function  and  constraints  and  then 
solving  the  resulting  Hamiltonian  boundary  value  problem.  In  contrast,  the  direct 
method  parameterizes  the  functions  using  function  approximation,  and  transcribes  the 
optimal  control  problem  to  a  nonlinear  programming  problem.  The  next  section 
examines  some  of  the  techniques  used  to  solve  dynamic  optimal  control  problems. 

2.1.1  Calculus  of  Variations 

The  Calculus  of  Variations  or  Method  of  Lagrange  multipliers  is  an  indirect 
method  in  which  the  dynamics  and  constraints  are  adjoined  to  the  performance  index 
with  a  time  varying  Lagrange  multiplier,  A  (i) .  The  new  cost  function  then  takes  the 
form: 

7  =  <i>(x(t0),t0,x(tf),tf)J  +  '  [ L(x(t),u(t),t)  +  AT  (t)  f  (x(t)  ,u(t)  ,t)-  x(t)jdt  (2.6) 
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A  stationary  solution  exists  when  dJ  =  0  for  arbitrary  du(t) .  In  other  words  the  cost  is 

either  at  a  minimum,  maximum  or  inflection  point,  and  the  formulation  of  the  cost 
function  determines  which  extreme  has  been  found.  From  calculus  of  variations  we  get 
the  Euler-Lagrange  equations  which  are  the  first-order  necessary  conditions  for  a 
stationary  solution  to  Equation  (2.1)  and  (2.6)  [11]. 

r)_  dHT  =  dLT  df  x 
c)x  c)x  dx 


(2.7) 


with  boundary  conditions 


(2.8) 


where 


H  (x(t)  ,u(t)  ,A(t)  ,t)  =  L(x(t) ,u(t) ,t)  +  AT  (t)  f  (x(t) ,u(t) ,t)  (2.9) 

and  a  stationary  solution  exists  if  the  optimality  criterion  is  satisfied: 

c)H 

—  =  0  ,t0<t<tf  (2.10) 

ou 

The  problem  reduces  to  a  Two-Point  Boundary-Value  Problem  (TPBVP)  where  the 
initial  state  and  final  Lagrange  multipliers  are  known.  The  optimal  solution  is  obtained 
using  a  'shooting'  algorithm  by  integrating  the  states  forward  in  time  and  then 
integrating  the  Lagrange  multipliers  backward  in  time.  The  resultant  Lagrange 
multipliers  are  called  costates.  The  costates  and  state  history  are  then  used  to  extract  the 
required  control  input.  This  method  has  a  "small  radii  of  convergence,"  and  the 
solution  may  become  unstable  and  diverge  while  either  integrating  the  states  forward  or 
integrating  the  Lagrange  multipliers  backwards  [26].  Lor  convergence,  this  method 


8 


requires  a  good  initial  guess  for  the  solution,  which  in  some  circumstances  is  difficult  to 
obtain. 


It  is  also  possible  to  employ  algorithms  that  look  at  the  gradient  of  the  dynamics 
and  the  performance  function.  Algorithms  like  Steepest  Descent  and  other  variations  on 
following  the  gradient  to  a  minimum  are  prone  to  only  finding  the  local  minimum.  In 
the  same  vein,  linearization  provides  a  means  of  determining  optimality  about  the  point 
of  expansion. 

2.1.2  Karush-Kuhn-Tucker  (KKT)  Conditions 

As  a  generalization  to  the  method  of  Lagrange  multipliers,  Karush,  Kuhn,  and 

Tucker  came  up  with  necessary  and  sufficient  conditions  for  optimal  solutions,  x°  (t ) 

where  the  superscript  'o'  denotes  optimal  [12, 15, 16]. 

The  conditions  are: 

1.  x°  (t) is  feasible 

2.  There  exist  A  ( t )  such  that: 

V7(jc0)  =  A(i)VC(x°) 

A(t)c(x°^  =  0 

A(t)  <  0 

The  KKT  conditions  provide  the  fundamental  tie  between  the  indirect  and  direct 
methods. 

2.1.3  Sequential  Quadratic  Programming  (SQP)  Methods 

Sequential  Quadratic  Programming  is  an  iterative  method  in  which  the  quadratic 
problem  is  solved  at  each  iteration.  Solving  the  quadratic  problem  at  each  iteration  will 
define  a  step  towards  the  next  search  direction  until  the  optimal  solution  is  found.  A 


(2.11) 
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quadratic  problem  consists  of  a  quadratic  model  of  the  cost  function  and  linearized 
constraints.  The  robustness  and  fast  convergence  of  SQP  methods  has  led  to  its 
incorporation  into  many  commercial  Nonlinear  Program  (NLP)  solvers  such  as  the 
industry  standard  Sparse  Nonlinear  Optimizer  (SNOPT)  from  Stanford  [17], 

2.1.3. 1  Nonlinear  Programming  (NLP) 

Nonlinear  programming  refers  to  a  specific  class  of  dynamic  optimization 
problems  in  which  the  objective  function  and/  or  constraints  are  nonlinear  and  has  the 
form: 


min  f(z) 
s.t.  g(z)<0 
h(z)  =  0 

zeZ 

where  /  is  the  objective/ cost  function  to  be  minimized,  g  is  the  inequality  constraints. 


(2.12) 


and  h  is  the  equality  constraints,  z  is  a  vector  of  parameters  that  exist  in  the  domain  of 
Z. 


Such  problems  often  require  the  equations  to  be  discretized  and  processed  in  a 
computer  algorithm.  The  name  'program'  in  nonlinear  programming  refers  to  its 
dependence  on  software.  Numerous  NLP  solvers  exist  such  as  SNOPT,  NPSOL,  and 
Interior  Point  Optimizer  (IPOPT).  Each  of  the  solvers  has  unique  features  that  make 
them  more  or  less  attractive  depending  on  the  problem  considered.  SNOPT  "employs  a 
sparse  SQP  algorithm  with  limited-memory  quasi-Newton  approximations  to  the 
Hessian  or  Lagrangian,"  whereas  NPSOL  employs  a  dense  SQP  algorithm  [4] .  Both 
SNOPT  and  NPSOL  are  effective  for  large-scale  nonlinear  problems  whose  functions 
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and  gradients  are  computationally  expensive  to  evaluate  [4] .  IPOPT  implements  an 
interior  point  method  for  solving  large-scale  nonlinear  optimization. 

2.2.4  Pseudospectral  Methods 

Pseudospectral  methods  are  a  direct  method  of  solving  the  optimal  control 
problem.  The  continuous  differential  equations  are  discretized  and  transcribed  to  a  NLP 
and  solved  by  one  of  many  available  NLP  solvers.  More  specifically,  it  is  a  form  of 
direct  collocation  where  the  state  and  its  derivative  are  both  determined  at  a  finite 
number  of  points.  The  problem  is  divided  into  segments  or  phases  and  the  state  is 
approximated  over  each  segment.  A  sample  schematic  of  segment  or  phase  linkages  is 
shown  in  Figure  2.  Phases  can  be  connected  in  any  order  and  may  have  multiple 
adjoining  phases.  The  approximation  is  differentiated  and  must  match  the  true  function 
derivative  at  a  finite  number  of  points  to  satisfy  the  collocation  conditions  and  ensure  a 
smooth  connection  between  the  phases  [26]. 


Figure  2.  Schematic  of  linkages  for  multiple-phase  optimal  control  problem  [25]. 
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The  states  are  approximated  using  one  of  three  collocation  schemes,  which  is  also 
the  primary  distinguishing  feature  between  the  various  pseudospectral  methods.  The 
three  collocation  schemes  are  Lobatto,  Radau,  and  Gauss.  Lobatto  collocation  schemes 
impose  the  collocation  condition  at  both  endpoints.  Examples  of  this  method  are 
Trapezoidal  (order  2)  and  Hermite-Simpson  (order  4)  quadratures.  Radau  collocation 
schemes  impose  the  collocation  condition  at  only  one  end  of  the  interval.  Gauss 
schemes  impose  the  collocation  conditions  strictly  interior  to  the  interval,  of  which  the 
simplest  Gauss  collocation  scheme  is  the  midpoint  rule  [9] . 

Each  of  the  segments  contains  a  finite  number  of  collocation  points  or  nodes. 

The  collocation  points  are  usually  the  roots  of  an  orthogonal  polynomial  although  this  is 
not  a  strict  requirement.  By  choosing  nodes  that  are  roots  of  an  orthogonal  polynomial, 
the  nodes  become  spaced  in  such  a  way  that  the  Runge  phenomenon  is  suppressed.  The 
roots  of  differing  degree  polynomials  are  unique  (except  at  zero  and  the  endpoints) 
thereby  ensuring  uniform  coverage  over  the  interval  of  interest.  The  nodes  and 
segments  (state  approximations)  are  then  tied  together  using  Lagrange  polynomials. 

The  theory  behind  pseudospectral  methods  is  well  documented  in  [22,  27]. 

There  are  a  few  software  packages  that  integrate  directly  into  MATLAB,  such  as 
DIDO  and  GPOCS,  which  use  pseudospectral  methods  to  solve  the  nonlinear  optimal 
control  problem.  DIDO,  written  by  I.  Michael  Ross,  uses  the  Legendre-Gauss-Lobatto 
pseudospectral  method  or  Lobatto  Pseudospectral  Method  (LPM)  [29],  DIDO  is  a 
software  package  that  was  initially  used  by  Maj.  Jorris  in  his  analysis  of  the 
maneuvering  re-entry  vehicle,  but  Maj.  Jorris  later  switched  over  to  Gauss 
Pseudospectral  Optimal  Control  Software  (GPOCS)  written  by  Anil  Rao  [25] .  GPOCS 
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uses  the  Legendre-Gauss  pseudospectral  method  or  Gauss  Pseudospectral  Method 
(GPM).  The  switch  from  DIDO  to  GPOCS  was  in  part  due  to  the  availability  of  costate 
and  Hamiltonian  information  output  by  GPOCS  that  directly  mapped  to  the  Lagrange 
multipliers.  This  was  a  key  feature  that  Maj.  Jorris  used  to  verify  the  optimality  of  the 
results  for  this  class  of  problems.  In  his  problem,  and  the  one  researched  herein,  the 
Hamiltonian  of  the  optimal  solution  is  equal  to  a  constant. 

H=- 1  (2.13) 

Figure  3  shows  the  equivalence  of  direct  and  indirect  methods  using  the  Gauss 
pseudospectral  discretization.  The  upper  left  box  begins  with  a  definition  of  the 
continuous-time  optimal  control  problem  such  as  that  presented  in  Section  2.1.  Using  an 
'indirect'  method  requires  formulating  the  first-order  optimality  conditions,  and  then 
solving  the  resulting  Hamiltonian  boundary-value  problem  (HBVP).  A  'direct'  method 
parameterizes  the  functions  using  function  approximation,  and  then  transcribes  the 
optimal  control  problem  to  a  nonlinear  programming  problem  (NLP)  to  solve  for  the 
parameters.  The  equivalence  of  the  two  solution  methods  are  then  shown  to  be  equal  by 
using  the  costate  mapping  theorem.  The  theory  behind  Gauss  Pseudospectral  methods 
can  be  found  in  [22,  25]  and  has  origins  in  Reddien's  paper  in  SIAM  Journal  on  Control 
and  Optimization  [26], 
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Figure  3.  Equivalence  of  indirect  and  direct  forms  using  the  Gauss  pseudospectral 

discretization  [25]. 


Much  of  this  research  builds  on  the  previous  work  of  Maj.  Jorris'  dissertation, 

[23] .  He  demonstrated  and  proved  that  a  direct  collocation  method  such  as  that  used  in 
pseudospectral  methods  could  be  used  to  find  optimal  solutions  for  the  re-entry 
problem  considered  herein.  The  pseudospectral  method  is  not  specific  to  a  particular 
application  such  as  re-entry  trajectory  generation,  rather  it  is  a  method  to  find  an 
optimal  solution. 

2.1.5  Classical  and  Modern  Control  Methods  (Inner-loop  Control) 

Although  not  addressed  by  this  research,  once  a  trajectory  is  obtained,  the  task  of 
following  the  optimal  trajectory  rests  on  the  inner-loop  controller.  The  inner-loop 
controller  may  include  a  human  being  or,  for  the  autonomous  system,  a  microprocessor 
alone.  The  optimal  trajectory  supplied  by  dynamic  optimization  methods  is  only  as 
good  as  the  inner-loop  control  tracking  the  reference  trajectory.  Inner-loop  control 
tracking  systems  include  models  of  the  dynamics  that  are  used  to  estimate  the  state  and 
provide  controls  to  drive  the  error  between  the  desired  trajectory  and  estimated 
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trajectory  to  zero.  Extensive  literature  and  examples  exist  for  linear  systems  to  track  an 
input;  however,  many  of  the  techniques  associated  with  linear  systems  cannot  be 
applied  to  the  nonlinear  re-entry  problem.  The  solution  is  to  transform  the  nonlinear 
equations  into  linear  equations  so  that  the  well  developed  linear  tools  can  be  used. 
Methods  for  transforming  the  problem  can  be  found  in  [30].  For  this  research,  the 
availability  of  suitable  inner-loop  control  is  assumed. 

2.2  Sensitivity  Methods 

"The  objective  of  sensitivity  analysis  is  to  quantify  the  effects  of  parameter 
variations  on  calculated  results"  [13].  This  research  seeks  to  characterize  how  changes  in 
waypoint  location  affect  the  terminal  state  (the  primary  mission).  Different  techniques 
and  theories  are  considered  for  applicability  to  solving  this  problem.  Included  in  the 
review  of  sensitivity  methods  are  perturbation,  sensitivity,  uncertainty  theory,  and  a 
direct  approach.  Particular  caution  is  needed  because  there  are  subtle  differences 
between  perturbation,  sensitivity,  and  uncertainty  analysis,  and  often  these  terms  gets 
mixed  in  the  literature. 

2.2.1  Perturbation  and  Sensitivity  Theory 

The  traditional  use  of  perturbation  and  sensitivity  theory  is  to  characterize  the 
stability  of  a  system  in  the  presence  of  perturbations  or  disturbances.  In  general  there 
are  three  basic  system  responses:  the  system  goes  unstable;  the  system  reacts  to  the 
disturbance,  but  remains  stable;  or  the  system  is  unaffected  by  the  disturbance.  A 
simple  example  is  the  pendulum,  which  exhibits  both  a  stable  and  unstable  mode.  The 
classic  pendulum  with  friction  when  subjected  to  a  disturbance  will  eventually  settle 
back  to  its  stable  equilibrium  point,  and  the  inverted  pendulum  will  move  away  from 
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the  equilibrium  point  due  to  the  presence  of  any  disturbances.  This  stability  theory, 
however,  while  noted,  is  not  the  primary  concern  of  this  research,  nor  what  is 
specifically  sought  after. 

2.2.2  Uncertainty  Analysis 

"The  objective  of  uncertainty  analysis  is  to  assess  the  effects  of  parameter 
uncertainties  on  the  uncertainties  in  calculated  results"  [13].  In  other  words,  uncertainty 
analysis  is  the  study  and  characterization  of  a  system  where  noise  has  been  added  to  the 
system.  This  is  an  important  aspect  in  inner-loop  control,  but  trajectory  generation  deals 
in  absolutes  precluding  this  as  an  option. 

2.2.3  Direct  Approach 

The  direct  approach  is  the  oldest  and  most  basic  tool  used  in  sensitivity  analysis. 
It  is  also  the  least  efficient  as  it  requires  solving  the  optimal  control  problem  repeatedly 
over  the  variations  or  perturbations  of  interest.  Due  to  the  inapplicability  of  the 
previous  methods,  this  method  serves  as  the  primary  tool  in  this  research.  The  direct 
approach  is  sufficient  for  mapping  the  sensitivity  of  an  optimal  control,  but  it  is  also  the 
most  computationally  intensive. 

2.3  3-D  Dynamics  Model  and  Fidelity 

The  previous  research,  from  which  this  investigation  is  derived,  used  a  3-D 
dynamics  model  with  the  following  assumptions  [23]: 

1.  Flat,  non-rotating  Earth. 

2.  Gravity  is  constant. 

3.  Flight-path  angle  is  small. 

4.  Drag  is  the  dominant  deceleration  term. 
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5.  Coefficient  of  Lift  (CL )  and  Coefficient  of  Drag  ( CD )  are  functions  of  angle-of- 
attack  only. 

Using  the  above  assumptions,  the  3-D  equations  of  motion  are: 


x  =  V  cos# 

(2.14) 

y  =  V  sin# 

(2.15) 

ll 

(2.16) 

t>  BV2e~Hh(l  +  cf) 

V~  2E* 

(2.17) 

f=BVe~^hc,  cos  C7  — —  +  U 
'  V 

(2.18) 

9  =  BVe~/3r°hcl  sin  a 

(2.19) 

where  x  and  y  are  lateral  displacements,  h  is  altitude,  V  is  velocity, /is  flight-path-angle, 
and  #  is  heading.  And  the  terms  on  the  right-hand  side  of  the  equal  sign:  B  is  a 
grouping  of  vehicle  parameters  and  initial  conditions,  (7  is  bank  angle,  C)  is  the  fraction 

of  CL  (coefficient  of  lift  that  results  in  the  maximum  lift-to-drag  ratio),  E  is  the 

maximum  lift-to-drag  ratio,  (3  is  the  atmospheric  scaling  height,  and  r0  is  a  reference 

radius.  The  derivation  and  further  explanation  of  the  terms  in  the  equations  can  be 
found  in  [23]. 

The  assumptions  above  are  suitable  for: 

1.  A  short  duration  trip  where  the  rotation  of  the  Earth  is  not  a  significant  factor. 

2.  A  short  range  trip  where  the  curvature  of  the  Earth  is  not  a  significant  factor. 

3.  Shallow  flight-path  angle  whereby  gravity  does  not  significantly  change. 

4.  A  vehicle  that  does  not  have  retro  firing  thrusters  or  similar  braking  forces. 
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The  assumptions  above  were  used  not  for  accuracy  in  modeling  reality,  but  for 
simplifying  the  dynamics  such  that  analytical  solutions  can  be  obtained  and  compared 
to  the  numerical  results.  Since  this  research  looks  to  increase  the  model  fidelity,  the 
assumptions  presented  above  will,  in  the  majority,  be  modified.  The  flat  Earth  will  be 
replaced  with  a  spherical  Earth.  Rotation  of  the  Earth  will  be  turned  back  'on.'  Gravity 
will  follow  an  inverse  square  law  as  a  function  of  radial  distance  from  the  center  of  the 
Earth.  Drag  will  remain  the  dominant  deceleration  term.  The  only  assumption  above 
that  will  be  simplified  further  is  that  the  coefficient  of  lift  and  drag  will  be  held  constant. 
In  general,  coefficients  of  lift  and  drag  are  a  function  of  angle-of-attack,  but  they  are  held 
constant  to  simplify  the  control  to  one  variable  (bank  angle). 

2.4  Highlighted  Vehicles 

Although  the  recent  GS  and  GPA  concepts  are  motivating  new  hypersonic 
research,  it  is  important  to  acknowledge  the  contributions  from  past  and  current 
contributions  of  hypersonic  vehicles.  Before  delving  into  vehicle  specifics,  some  classic 
definitions  follow  in  the  next  section. 

2.4.1  Vehicle  Descriptors 

The  current  buzzwords  or  acronyms  describing  hypersonic  vehicles  are  CAV, 
HCV,  and  RLV.  The  terms  change  over  time,  but  the  underlying  concepts  remain  the 
same.  The  definition  of  each  term  is  now  presented. 

2.4.1. 1  Common  Aero  Vehicle  (CAV) 

The  Common  Aero  Vehicle  (CAV)  is  a  concept  which  describes  a  space 
re-entry  aeroshell  launched  into  space  on  a  suitable  vehicle,  which  then 
survives  atmospheric  re-entry,  reduces  its  speed  to  low  Mach  numbers  or 
even  sub-Mach,  and  dispenses  a  cargo  payload  or  weapon  in  the  Earth's 
atmosphere.  The  conceptual  CAV  might  be  propelled  into  space  by  any 
of  a  number  of  present  or  future  launch  platforms  and  dispense  a  wide 
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variety  of  cargoes,  payloads  or  weapons.  Development  of  the  CAV 
capability  will  satisfy  future  requirements  enunciated  in  numerous 
national  visions,  future  studies,  and  military  plans  and  will  ultimately  be 
necessary  to  fully  realize  the  opportunities  inherent  in  operating  from, 
through,  and  in  space  and  give  true  meaning  to  the  phrases  global  reach, 
global  power  projection,  and  global  engagement.  [28] 

2. 4.2. 2  Hypersonic  Cruise  Vehicle  (HCV) 

HCV  describes  a  generic  class  of  vehicles  capable  of  maintaining  hypersonic 
speeds  under  thrust  for  sustained  periods  of  time.  At  the  moment,  HCVs  are 
conceptual,  but  include  several  projects  such  as  MA3A,  DARPA  sponsored  FALCON, 
Black-swift,  and  X-51.  The  US  military  is  seeking  an  unmanned  HCV  capable  of 
delivering  12,000  lbs  of  payload  that  can  reach  any  target  within  a  range  of  9,000  nautical 
miles  in  less  than  two  hours  [2,  23]. 

An  extension  of  this  research  directly  applies  to  HCVs,  specifically  the  addition 
of  thrust  to  the  equations  of  motion.  The  HCV  that  burns  out  of  thrust  is  subject  to  the 
same  dynamics  as  the  CAV. 

2.4.1.3  Reusable  Launch  Vehicle  (RLV) 

The  term  RLV  comes  from  NASA's  Reusable  Launch  Vehicle  program.  The 

program  is  now  into  its  Second  Generation  RLV  program,  and  is  also  known  as  the 

Space  Launch  Initiative.  The  goals  of  the  program  are: 

NASA's  Reusable  Launch  Vehicle  (RLV)  program  has  dual  objectives:  to 
demonstrate  technologies  leading  to  a  new  generation  of  space  boosters 
capable  of  delivering  payloads  at  significantly  lower  cost,  and  to  provide 
a  technology  base  for  development  of  advanced  commercial  launch 
systems  that  will  make  U.S.  aerospace  manufacturers  more  competitive  in 
the  global  market.  [5] 

An  RLV  would  then  only  require  refueling  for  subsequent  returns  to  space. 
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RLVs  are  considered  important  beneficiaries  of  this  research.  In  general,  an  RLV 
will  either  need  to  fly  back  to  its  home  designation  or  maneuver  to  another  facility 
downrange.  The  decision  to  turn  around  and  fly  home  or  continue  on  to  another 
landing  site  is  contingent  upon  the  time  of  payload  separation  and  whether  the  vehicle 
has  enough  energy  and  control  authority  to  make  the  return  trip. 

2.4.2  X-15 


In  1954  the  X-15  became  the  first  suborbital  airplane.  The  X  designation  denotes 
the  plane  as  an  experimental  vehicle,  but  also,  and  more  importantly  as  a  high  altitude 
research  vehicle.  It  significantly  "contributed  valuable  research  information  in  the 
supersonic  and  hypersonic  speed  regime  up  to  the  fringes  of  space"  and  expanded  our 
understanding  and  modeling  of  the  atmosphere  [24],  The  X-15  is  not  capable  of  taking 
off  from  the  ground  on  its  own;  therefore  it  is  attached  to  a  B-52  which  will  release  it 
midair.  Figure  4  shows  the  X-15  separating  from  the  B-52. 


Figure  4.  This  photo,  taken  from  a  B-52,  pictures  the  X-15  immediately  after  launch  with 

an  F-104  flying  chase  [7], 
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2.4.3  Dynamic  Soarer  (Dyna-Soar) 

Perhaps  the  first  CAV-like  design,  the  Dyna-Soar,  also  known  as  the  X-20,  was  a 
proposed  spaceplane.  It  was  designed  to  be  a  multi-purpose  vehicle  that  would  be 
boosted  to  orbit  and  then  glide  to  its  destination.  In  addition  to  orbital  bomber  and 
reconnaissance  missions,  designers  also  added  satellite  maintenance  and  sabotage 
missions  to  its  list  of  capabilities.  Ultimately  the  project  was  cancelled  in  1963  due  to 
budget  overruns  and  technical  hurdles.  Although  this  project  never  got  off  the  ground, 
the  research  and  design  concepts  made  its  way  into  other  projects  such  as  the  Space 
Transport  System  (STS),  and  many  of  the  planned  capabilities  are  now  incorporated  in 
the  Air  Force's  current  Global  Reach  concepts  [21], 


Figure  5.  Boeing's  mockup  of  the  X-20  [3]. 


2.4.4  Space  Transport  System  (STS) 

STS,  better  known  as  NASA's  Space  Shuttle,  could  be  considered  an  RLV  and  a 
CAV.  The  STS  includes  the  orbiter,  solid  rocket  boosters  and  external  fuel  storage  tank. 
The  solid  rocket  boosters  and  external  tank  are  jettisoned  after  launch,  but  the  boosters 
parachute  down  for  recovery  and  reuse  (unlike  the  external  tank).  Additionally,  the 
Shuttle's  return  to  Earth  begins  as  a  hypersonic  unpowered  glide  which  is  similar  in 
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description  to  the  CAV,  but  unlike  the  CAV  the  Shuttle  is  not  dispensing  a  cargo  or 
payload  while  in  the  atmosphere. 


Figure  6.  Discovery  shuttle  and  STS  120  crew  on  launch  to  the  International  Space 

Station  [6], 

2.4.5  Flexible  Aerospace  System  Solution  for  Transformation  (FASST) 

"The  Boeing  FASST  concept  is  a  horizontal  takeoff,  horizontal  landing  (HTHL), 
two-stage-to-orbit  (TSTO)  reusable  vehicle  with  a  turbine-based  propulsion  system  on 
the  first  stage  and  a  rocket-based  combined  cycle  (RBCC)  powered  second  stage"  [10]. 

FASST  is  a  response  to  NASA's  Integrated  Space  Transportation  Plan/Program 
(ISTP)  which  seeks  to  expand  the  civil  and  commercial  reach  into  space  for  the  coming 
decades.  The  plan  details  the  development  of  RLVs  for  use  in  delivering  cargo  to  Low 
Earth  Orbit  (LEO).  The  turbine  stage  of  the  vehicle  is  set  for  an  upper  limit  of  Mach  4 

[19]. 
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2.4.6  Force  Application  Launch  from  Continental  United  States  (FALCON) 

FALCON  is  currently  a  DARPA  sponsored  project  whose  "objectives  are  to 
develop  and  demonstrate  hypersonic  technologies  that  will  enable  prompt  global  reach 
missions"  [2],  The  project  expands  upon  the  path  of  developing  a  HCV  and  has  already 
spawned  several  sub-projects  labeled  as  Hypersonic  Technology  Vehicle(s)  (HTV).  This 
test  bed  of  vehicles  will  be  used  to  validate  affordable  and  responsive  space  lift 
capabilities. 

2.5  Closing  Remarks 

This  chapter  presented  the  dynamic  optimization  problem  and  briefly  discussed 
techniques  used  to  solve  the  sensitivity  problem.  The  verification  of  the  optimality  of 
the  results  obtained  from  the  pseudospectral  solver  is  expounded  upon  in  Maj.  Jorris' 
dissertation  [23],  Additionally,  several  real  world  vehicles,  to  which  this  research  may 
potentially  apply,  were  introduced. 

During  the  course  of  this  research,  the  software  written  by  Anil  Rao,  and  used 
herein,  changed  names  from  GPOCS  to  OPENPOCS  to  GPOP  to  PSCOL.  Each  name 
change  signifies  either  a  change  in  code  structure,  licensing,  or  product  affiliation,  but 
the  software  maintained  the  underlying  GPM.  In  its  current  release,  PSCOL  can  also  do 
both  the  LPM  and  Radau  Pseudospectral  Method  (RPM).  This  research  uses 
OPENPOCS  for  legacy  reasons,  although  any  pseudospectral  method  may  be 
considered. 
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III.  Problem  Definitions  and  Assumptions 

3.1  Generic  Problem  Statement 

The  overall  objective  of  this  research  is  to  obtain  the  optimal  trajectory  for  a 
given  vehicle  to  reach  its  target  in  a  minimum  amount  of  time  and  to  study  the 
sensitivity  of  the  terminal  state  due  to  variations  in  waypoint  location(s).  The  optimal 
solution  to  the  two-point  boundary-value  problem  will  serve  as  the  baseline  or  nominal 
solution  to  which  comparisons  will  be  made. 

Additionally,  a  verification  of  the  need  for  a  higher  fidelity  model  is  sought. 
Without  argument,  people  accept  that  the  Earth  is  spherical  (the  Earth  is  actually  an 
oblate  spheroid),  and  a  flat  Earth  model  is  insufficient  for  any  orbit  problem  including 
re-entry;  therefore,  without  further  consideration,  the  flat  Earth  model  is  rejected  in  lieu 
of  the  spherical  Earth  model.  The  verification  then  rests  on  the  need  to  include  Earth's 
rotation  in  the  model. 

3.2  Mission  Assnmptions 

The  assumptions  are  broken  into  three  classes:  mission  components,  vehicle 
components,  and  environmental  components.  Mission  components  include  target  and 
waypoints.  Vehicle  components  include  equations  of  motion  and  control  limitations. 
3.2.1  Mission  Components 

Similar  mission  assumptions  are  taken  from  [23]: 

1.  The  waypoints  are  specified  in  the  desired  sequence. 

2.  Waypoint  passage  must  be  through  a  radial  originating  from  the  center  of  the 
Earth. 

3.  Altitude  for  waypoint  passage  is  not  specified. 
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4.  Inner-loop  control  is  available.  Only  the  outer-loop  or  trajectory  generation  is 
addressed. 

5.  Target  coordinates  and  final  altitude  is  specified. 

3.2.2  Vehicle  Components 

1.  Vehicle  is  modeled  as  a  point  mass. 

2.  Non-thrusting  vehicle. 

3.  Bank  angle  is  the  only  control  parameter. 

4.  Coefficients  of  lift  and  drag  are  held  constant. 

3.2.3  Environmental  Components 

1.  Spherical  Earth  model. 

2.  The  atmosphere  is  modeled  as  a  simple  exponential. 

3.  Rotation  of  the  Earth  is  turned  'off'  initially. 

3. 2.3.1  3-D  Non-dimensional  Equations  of  Motion 

In  numerical  dynamic  optimization,  it  is  desirable  to  scale  and /  or  non- 

dimensionalize  the  equations  of  motion  such  that  all  the  terms  involved  are  near  the 

same  order  of  magnitude.  Vinh  and  Braces'  non-dimensional  "universal  equations"  of 

motion  are  problematic  because  they  are  not  scaled  appropriately  for  dynamic 

optimization  solvers.  Instead,  the  dimensional  3D  equations  of  motion  from  [20]  are 

scaled  and  non-dimensionalized.  The  following  are  the  3D  non-dimensional  equations 

of  motion.  The  derivation  of  these  equations  is  provided  in  Appendix  A. 

dr  R  ,  . 

—  =  V  sin  y 
dr 

dRV'  sin  y  ,  ,2  ,,  ,  .  ...  s 

- = - — — b  r  (Of  cos^(cos0sin  y-suupsvny/ cosy)- D 

dr  r 


(3-1) 

(3.2) 
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eld'  _  f  RV'  |  cos  ycos  y 
dr  y  r  J  cosp 


d(j>  f  V^| 

- =  — —  cosy  sin  ^ 

dr  y  r  J 


(3.3) 

(3.4) 


dy  L'  cos  y  RV'  , 

—  =  — 008  0---^^  +  — COS7+2«e  cos^cos^  +  ... 


. . .  +  -^—7  cdk ’l  cos  <j)  ( cos  <j)  cos  y  +  sin  (j)  sin  y  sin  y) 


dy'  L' sine  RV'  ,  _  ,  ,  .  ,  . 

— —  =  - - —  cos  ^cos  y  tan  ep  +  2 (Oe  ( sin  y  cos  ep  tan  y -  sin  ep) . . . 

dr  V  cos  y  r 

ret)#  sin (Zi cos ^ cos y 

V  cos  y 

where  r  is  the  radius  from  the  center  of  the  Earth,  RV  is  the  relative  velocity,  0  is  the 
longitude,  <p  is  the  latitude,  y  is  the  flight-path  angle,  and  y  is  the  heading  angle.  The 
terms  (not  already  defined)  on  the  right-hand  side  of  the  equal  sign:  L  is  the  force  due 
to  lift,  D  is  the  force  due  to  drag,  and  <y0  is  the  rotation  rate  of  the  Earth.  The  prime 


denotes  the  values  as  non-dimensional. 


3.2.4  Mission  Profiles  and  Case  Definitions 

Several  mission  profiles  are  considered  in  which  the  RV  will  attempt  to  fly  a 
minimum  time  trajectory  to  the  target.  The  cases  are  designed  to  investigate  the  effect  of 
Earth's  rotation  on  generated  trajectories,  and  the  last  case  examined  is  an  example  of 
how  a  mission  strategist  might  use  this  tool  to  select  additional  waypoints.  Each  of  the 
cases  share  the  same  initial  and  final  target  coordinates,  but  otherwise  exhibit  a 
progression  of  changes.  The  first  mission  profile  (case  1)  is  used  to  validate  the 
equations  and  solution  technique,  and  it  begins  with  a  scenario  in  which  the  rotation  of 
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the  Earth  is  turned  'off.'  This  scenario  is  then  contrasted  against  the  second  mission 
profile  (case  2)  where  the  rotation  is  turned  'on.'  Although  Table  1  and  Table  2  show  the 
two  cases  with  identical  initial  conditions,  in  reality  the  initial  conditions  are  different 
due  to  inertial  and  relative  measurements  of  the  velocity  and  flight-path  angle.  The 
superscript  'R'  in  the  tables  is  used  to  denote  a  relative  value.  The  third  mission  profile 
(case  3)  is  an  extension  of  case  2  where  an  intermediate  waypoint  has  been  defined.  In 
all  cases,  the  optimal  trajectory  is  computed  between  the  initial,  intermediate  (if 
defined),  and  final  target.  The  computed  optimal  trajectory  defines  the  baseline 
trajectory  from  which  the  sensitivity  analysis  stems.  The  sensitivity  analysis  is  explained 
in  Section  3.4.  Specific  values  for  each  case  are  presented  in  Table  1  to  Table  3 
respectively,  with  the  initial,  intermediate  (if  defined),  and  final  state  values  defined. 


Table  1.  Case  1  initial  and  final  conditions. 


Initial  Conditions 

Final  Conditions 

r0  =  6500 Am 

rf  =  645 6 Am 

V0  =7. 3152 

Vf=ND 

r0=- 1-5° 

Yj  =  ND 

3 

ll 

o 

Yf  =  ND 

e0  =  -80° 

ef  =ii.32o 

CO 

II 

</>f  =3.52° 

£70  =  ND 

a  f  =  ND 

<y0  =  0-^-,  m  =  lOOOAg  ,  S  =  80m2,  ND=not  defined 
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Table  2.  Case  2  initial  and  final  conditions. 


Initial  Conditions 

Final  Conditions 

r0  =  6500 km 

rf  =  6456 km 

\  =7.3152^ 

RVf  =  ND 

Rr0 =-i-5° 

Ryf  =  ND 

-3 

ii 

o 

y/f  =  ND 

II 

1 

OO 

o 

9f  =11.32° 

</>0  =  34” 

(/>f  =3.52° 

<Tq  =  ND 

<jf=ND 

<y0  =  7.2722x10  5  ^ ,  m  =  500 kg  ,  S  =  80m2,  ND=not  defined 

Table  3.  Case  3  initial,  intermediate,  and  final  conditions. 


Initial  Conditions 

Intermediate  Waypoint 

Final  Conditions 

r0  =  6500 km 
rV0  =  7. 3152 

RVo  =  _l-5° 

¥ o  =  °° 

90  =  -80° 

</>0  =  34° 

<J0  =  AT) 

rx  =  ND 

Ry1  =  nd 

Ry  =  ND 
y/x  =  ND 

9X  =  -22.07° 

$  =21.59° 
f tx=ND 

rf  =  645 6&m 
'V,  =  AT> 

V/  =  ND 
y/f  =  M) 
ef  =  ii.32o 
</>f  =  3.52° 

<jf=ND 

<y0  =  7.2722x10  5  ^ ,  m  =  500kg  ,  S  =  80m2,  ND=not  defined 

3.2.5  Vehicle  Characteristics 

The  vehicle  characteristics  chosen  are  not  specific  to  any  real  world  vehicle,  such 
as  presented  in  Chapter  II.  Instead  they  represent  "ballpark"  figures  that  describe  a 
small-to-medium  size  vehicle  (case  1),  and  a  similar  size  vehicle  with  half  the  weight 
(cases  2  and  3).  Since  the  initial  state  vector  is  held  constant,  a  difference  in  initial  kinetic 
energy  is  introduced  when  the  rotation  of  the  Earth  is  turned  on.  The  vehicle  weight  is 
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cut  in  half  to  reduce  the  initial  inertial  kinetic  energy  and  prevent  the  vehicle  from 
bouncing  off  the  Earth's  atmosphere  when  the  rotation  of  the  Earth  is  turned  'on.'  The 
vehicle  characteristics  are  simply  described  by  the  following  parameters: 

S  =  80 m2  -  wetted  surface  area 


m  ■ 


:  1000%  (500% )  -  mass  of  RV 


CL=  2-  coefficient  of  lift 
CD  =  1  -  coefficient  of  drag 

The  vehicle  lift  and  drag  coefficients  are  simplified  to  be  constants.  This  corresponds  to 
setting  the  angle-of-attack  to  a  constant,  since  the  lift  and  drag  are  primarily  a  function 
of  angle-of-attack.  This  is  consistent  with  the  defined  equations  of  motion  in  Section 
3. 2.3.1  which  only  include  bank-angle  as  a  control.  The  RV  is  given  a  minimum  and 
maximum  allowable  bank-angle.  Such  limitations  are  common  on  vehicles  that  are 
shielded  more  heavily  on  one  side.  Additionally  the  bank-angle  rate  is  constrained 
reflecting  that  the  vehicle  cannot  instantaneously  change  its  bank-angle. 

<7  .  =-60° 

min 

a  =60° 

max 

^min  =  -60  %ec 
^max=60  %ec 

3.3  Optimal  Control  Problem 

As  stated  earlier,  the  objective  is  to  minimize  the  time  of  flight;  therefore,  the 
dynamic  optimization  problem  is: 

Minimize 

J  =  7f 


(3.7) 


(3.8) 
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subject  to  the  non-dimensional  equations  of  motion,  (3.1)-(3.6)  as  well  as  the  dynamics, 
representing  the  relationship  between  the  bank-angle  rate  and  the  control  input  u  (t)  . 


LI  KJ  /  \ 

—  =  u(t)  (3.9) 

dr 

Adding  the  bank-angle  rate,  (3.7),  as  a  state  allows  for  limits  on  how  fast  the  vehicle  can 
roll.  The  TPBVP  is  subject  to  the  event  constraints  or  boundary  conditions,  specifically: 
initial  altitude,  velocity,  longitude,  latitude,  flight-path-angle,  and  heading,  as  well  as 
final  altitude,  longitude,  and  latitude.  These  are  given  as: 


r'(To)~ro 

RV'(r0)-V0 


$  (*■„) 

f(T0)-r0 

V'{?o)-Vn 

r'(Tf)-'r 


=  0 


(j)  [rf )  <pf 


(3.10) 


Additionally,  the  problem  is  subject  to  the  inequality  constraints  on  the  state  (bank- 
angle)  and  control  (bank-angle  rate). 


C(x(t),u(t),t) 


(7  • 

min  V  / 

g(t)-g 
&  —u(t) 
u(t)-g 

V  /  max 


(3.11) 


3.4  Methodology 

This  section  describes  the  systematic  selection  of  waypoints  used  to  map  the 
sensitivity  of  the  terminal  state.  Preliminary  experimentation  and  testing  helped  define 
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the  flowchart  in  Figure  7.  First,  a  single  phase  optimal  trajectory  is  generated  using 
(3.10)  and  (3.11)  but  without  any  interior  point  constraints  (waypoints).  This  trajectory 
serves  as  a  baseline  or  nominal  trajectory  for  future  computations.  Using  Bellman's 
principle  of  optimality,  the  single  phase  solution  is  divided  into  a  two-phase  solution. 
The  connection  point  of  the  two  phases  is  a  waypoint,  or  fixed  latitude  and  longitude  the 
RV  must  pass  through.  Splitting  the  trajectory  into  two  phases  now  provides  a 
mechanism  to  precisely  define  and  relocate  waypoints.  The  flow  chart  outlines  the 
progression  or  variation  of  the  intermediate  waypoints  in  a  systematic  fashion. 

The  flow  chart  shown  begins  with  user  supplied  initial  and  final  conditions. 
Specifically  altitude,  latitude,  and  longitude  are  needed  for  both  initial  and  final 
conditions,  and  in  addition  to  the  previously  mentioned  states,  velocity,  flight-path 
angle,  and  heading  angle  are  required  for  the  initial  condition.  The  remaining  sequence 
of  events  is  hands-off  for  the  user.  The  continuous  time  problem  is  discretized  and 
evaluated  at  n  LG  nodes  to  solve  the  baseline  optimal  trajectory  from  the  initial  to  final 
target.  Each  node  of  the  baseline  optimal  solution  then  becomes  a  waypoint,  and  at  this 
point  the  single  phase  trajectory  is  divided  into  many  two-phase  trajectories.  New 
optimal  control  problems  are  created  by  incrementally  shifting  the  waypoint  north  by 
'dphi'  degrees.  The  optimal  trajectory  is  found  through  each  new  waypoint  position, 
and  the  process  of  incrementing  by  'dphi'  continues  until  a  solution  cannot  be  found, 
'dphi'  was  chosen  to  be  a  relatively  small  value,  i.e.  'dphi'  =  0.01  (-0.6  nmi),  to  aid  in  the 
convergence  rate  of  the  GPM  solver  which  uses  the  nearest  computed  trajectory  as  an 
initial  guess  for  the  new  waypoint.  Larger  values  of  'dphi'  may  be  used,  but  it  increases 
the  likelihood  the  GPM  solver  will  not  converge  since  the  initial  guess  is  comparatively 
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Figure  7.  Basic  algorithm  used  to  map  sensitivity. 

farther  from  the  optimal  solution.  The  exact  specifics  of  what  constitutes  an 
unobtainable  solution  will  be  explained  later.  Next,  the  same  procedure  is  repeated. 
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incrementally  shifting  the  waypoints  south.  Note  that  the  solution  is  optimal  at  the 
nodes  and  it  is  assumed  that  the  solution  is  near  optimal  elsewhere.  Splines  are  used  to 
interpolate  between  the  nodes  and  provide  a  more  accurate  solution.  Tabulating  the 
terminal  state  values  associated  with  each  waypoint  then  describes  how  the  terminal 
state  changes  with  waypoint  location,  providing  the  desired  'mission  objective' 
sensitivity  to  the  waypoint. 

Figure  8  visually  conceptualizes  the  process  in  the  flowchart.  First,  the  baseline 
optimal  solution  between  the  initial  and  final  target  is  computed  shown  by  the  solid  line. 
Then,  starting  from  a  location  on  the  baseline  optimal  trajectory,  waypoints  are 
incrementally  stepped  along  a  longitude  in  both  directions  in  a  uniform  increment.  At 
each  new  waypoint,  the  optimal  trajectory  through  the  new  waypoint  is  calculated 
shown  in  dashed  lines.  Once  a  trajectory  is  considered  unobtainable,  the  process  would 
then  repeat  by  starting  over  at  a  new  location  along  the  baseline  trajectory. 

There  are  several  automated  'stop'  conditions  in  the  software  code  that  are  used 
to  manage  the  iteration  control  and  convergence  criteria  the  program  uses  to  obtain  a 
solution.  A  'stop'  condition  is  a  message  to  the  program  to  stop  stepping  the  waypoint 
north  or  south,  and  move  on  to  the  next  waypoint  location  on  the  baseline  solution. 

Since  this  research  is  primarily  interested  in  the  immediate  surrounding  area  of  the 
baseline  solution,  there  is  a  limit  on  the  number  of  incremental  steps  a  waypoint  may  be 
moved  north  or  south.  A  'stop'  condition  may  also  be  triggered  if  the  underlying 
optimization  solver  returns  back  status  codes  corresponding  to  infeasible  solutions  or 
other  unexpected  termination  status  codes.  Before  the  'stop'  condition  is  triggered,  the 
program  increases  the  number  of  nodes  in  each  phase  and  attempts  to  re-solve  the 
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problem.  Once  the  number  of  nodes  in  each  phase  has  reached  a  limit,  then  the  'stop' 


condition  is  triggered.  The  justification  for  increasing  the  nodes  is  presented  in  Section 
4.1. 


Figure  8.  Waypoint  selection  process  w.r.t.  optimal  trajectory. 


This  methodology  meets  the  objectives  of  the  research  by  solving  the  baseline 
optimal  control  problem,  and  then  sampling  the  neighboring  area  for  fixed  waypoints 
and  re-solving  the  optimal  trajectory  to  obtain  sensitivity  data.  The  expected  results  are 
a  range  of  feasible  solutions  that  fall  in  an  envelope  surrounding  the  baseline  trajectory. 
Specifically,  the  time-to-target  will  increase,  and  the  terminal  velocity  will  decrease  as 
the  RV  maneuvers  to  waypoints  farther  off  the  baseline  trajectory.  The  results  obtained 
using  this  approach  is  discussed  in  Chapter  IV. 
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IV.  Analysis  and  Results 


In  this  chapter,  the  analyses  of  the  three  cases  defined  in  Chapter  III  are  used  to 
present  a  unifying  trend  among  all  re-entry  trajectories.  Specifically,  cases  1  and  2  are 
similar  in  problem  setup  and  vehicle  dynamics,  but  they  differ  in  vehicle  mass  and  Earth 
rotation  rate.  The  last  case  is  more  of  a  demonstration  of  the  further  capabilities  of  this 
solution  technique.  The  numerical  result  obtained  from  each  of  the  cases  fulfills  the 
objective  of  identifying  the  sensitivity  of  the  terminal  state  to  variations  in  waypoint 
location. 

4.1  Node  Analysis 

The  amount  of  time  to  compute  an  optimal  solution  in  general  increases  as  more 
nodes  are  used  in  the  solution.  Recall  a  node  represents  a  discrete  time  in  the 
continuous  time  dynamics.  The  added  benefit  of  using  more  nodes  is  obtaining  a  more 
accurate  solution;  therefore,  the  motivation  for  this  portion  of  the  analysis  is  to  identify  a 
minimum  number  of  nodes  that  should  be  used  in  calculating  solutions,  balancing 
computing  time  and  ensuring  accuracy  in  the  solution. 

To  establish  a  true  baseline,  an  initial  state  vector  was  chosen  and  integrated 
forward  for  approximately  2000  seconds  using  MATLAB's  ode45  integrator.  The  control 
variable  (bank-angle)  was  fixed  at  0°  for  the  entire  trajectory.  Now  using  the  altitude, 
latitude,  and  longitude  at  2000  seconds  as  the  final  target  vector,  the  initial  and  final 
vector  are  fed  into  the  pseudospectral  solver  for  processing.  The  state  vector  and  vehicle 
parameters  supplied  to  the  program  are  listed  in  Table  4. 

Despite  the  pseudospectral  solver  returning  successful  messages  of  convergence 
to  an  optimal  solution,  the  pseudospectral  solver  generated  differing  trajectories  based 
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on  the  number  of  nodes  used.  Fortunately,  the  results  increased  in  accuracy  as  the 
number  of  nodes  increased,  i.e.  more  closely  matched  the  baseline.  To  elucidate, 
solutions  for  twenty-one  nodes  and  sixty-one  nodes  are  investigated. 

Table  5  shows  the  resultant  terminal  state  vector  obtained  from  using  an  ode45 
solver,  and  a  pseudospectral  solver  with  twenty-one  and  sixty-one  nodes  respectively. 
Notably  the  twenty-one  node  solution  differed  significantly  from  the  baseline  (ode45)  in 
both  terminal  state  conditions  and  the  amount  of  time  required  to  reach  the  target.  In 
contrast,  the  sixty-one  node  solution  nearly  matched  the  true  solution  with  exception  to 
the  final  bank  angle. 

Figure  9  shows  the  altitude  of  the  respective  solutions  with  respect  to  time.  A 
cubic  polynomial  is  fit  through  the  nodes,  represented  by  'x's.  Note,  the  blue  and  green 
lines  corresponding  to  the  sixty-one  node  and  ode45  solutions  are  indistinguishable 
from  each  other.  Graphically  it  is  quite  apparent  that  using  twenty-one  nodes  will  not 
suffice.  Towards  the  middle  of  the  trajectory  where  the  nodes  are  spread  the  farthest, 
the  solution  using  twenty-one  nodes  is  unable  to  capture  the  true  dynamics.  This  is  a 
consequence  of  the  sampling  rate  and  is  remedied  by  increasing  the  sampling  rate  or 
increasing  the  number  of  nodes  used. 
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Table  4.  Initial  and  final  state  vectors  supplied  to  the  optimal  solver. 


Initial  Conditions 

Final  Conditions 

r0  =  6500 km 

rf  =  6456A/77 

V0  =73152 if 

V f  ~  ND 

7o  =  -1-5° 

rf  =  ND 

ii 

o 

y/f  =  ND 

II 

1 

oo 

q 

ef  =n.32o 

</>o  =  34° 

(j)f  =  3.52° 

£70  =  ND 

o f  =  ND 

coffi  =  7.277x10  5  — ,  ND=not  defined 

©  sec  ' 

Table  5.  Three  terminal  state  vectors  are  shown  that  are  propagated  from  the  same 
initial  state  vector  (see  Table  4)  using  MATLAB's  ode45  and  a  pseudospectral  solver 
with  21  and  61  nodes  respectively. 


MATLAB's  ode45 

21  Nodes 

61  Nodes 

Final  Conditions 

rf  =  6456km 

rf  =  6456177? 

rf  —  6456 km 

Vf  =  1.592 

Vf  =1.702 

Vf  =1.592 

yf  =  -1.04° 

yf  =-17.32° 

yf  =-1.02° 

y/f  =  -31.04° 

y/f  =-28.51° 

y/f  =-31.17° 

Qf  =11.32° 

6f  =11.32° 

ef  =11.32° 

(j)f  =  3.52° 

<pf  =  3.52° 

</>f  =3.52° 

ii 

o 

(jf  =  60° 

af  =19.86° 

Final  Time 

2001  seconds 

2016  seconds 

2001  seconds 
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altitude  [km] 


Re-entry  Trajectory 


Figure  9.  A  re-entry  trajectory  solution  generated  using  21  /  61  nodes  in  OPENPOCS,  and 
a  Runge-Kutta  4-5  integrator. 
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Further,  Figure  10  shows  the  Hamiltonian  for  each  of  the  solutions.  Stated  earlier 
the  Hamiltonian  of  the  true  optimal  solution  is  equal  to  a  constant,  H  =  —  1 .  Clearly  the 
twenty-one  node  solution  is  not  optimal,  and  the  sixty-one  node  solution  is  nearly 
optimal.  Although  finding  an  optimal  minimum  number  of  nodes  for  all  trajectories  is 
beyond  the  scope  of  this  research,  it  is  clear  that  more  nodes  are  better.  For  the 
sensitivity  study  herein,  the  sixty-one  node  solution  will  be  considered  adequate; 
therefore,  for  further  calculations,  a  minimum  of  sixty-one  nodes  is  used  for  the  global 
solution. 


39 


4.2  Rotation  Effects 

The  inclusion  of  the  rotation  rate  of  the  Earth  is  a  significant  parameter  in 
modeling  reality.  Figure  11  shows  the  range  of  a  re-entry  vehicle  in  constant  zero- 
degree  bank-angle  flight  independent  of  initial  heading  angle  for  the  Earth  rotation  both 
'off'  and  'on.'  Note  that  the  figure  shown  is  a  flat  representation  of  the  spherical  Earth, 
e.g.  a  Mercator  map  projection.  To  ensure  a  true  comparison,  the  inertial  velocity 
(kinetic  energy)  and  flight-path  angle  of  the  RV  are  held  constant.  With  the  rotation  of 
the  Earth  'off',  the  RV  reaches  targets  in  all  directions  at  an  equal  distance  and  reaches 
the  ground  at  the  same  time  regardless  of  heading.  For  the  rotating  Earth,  the  relative 
velocity  and  flight-path  angle  will  vary  as  the  heading  angle  changes.  Consequently,  the 
solid  line,  depicting  Earth's  rotation  'on',  shows  that  the  RV  travels  farther  in  a  westerly 
direction  and  shorter  in  an  easterly  direction.  Note  that  the  RV  will  also  reach  the 
ground  at  different  times.  The  footprint  is  not  merely  shifted  to  the  west,  but  it  is  also 
slightly  skewed  due  to  the  fact  that  the  Earth  below  the  RV  moves  at  a  different  rate  as  a 
function  of  latitude.  It  should  be  noted  that  the  impact  footprint  size  is  a  function  of  the 
vehicle's  lift-to-drag  ratio  (L/D).  This  is  because  L/D  is  a  measure  of  aerodynamic 
efficiency,  and  in  general,  aerodynamic  efficiency  and  range  are  roughly  directly 
proportional.  Thus,  Figure  11  is  only  applicable  for  the  vehicle  defined  in  Table  2,  and 
would  need  to  be  recomputed  for  other  vehicle  configurations.  In  conclusion,  final 
targets  that  were  obtainable  with  the  stationary  Earth  model  may  now  be  infeasible  with 
the  rotation  turned  'on.' 

Figure  11  was  created  by  choosing  an  initial  inertial  state  vector,  and  then 
propagating  the  state  forward  in  time  until  the  RV  hit  the  ground  (altitude  equals  zero). 
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For  the  rotating  Earth,  the  initial  inertial  velocity  and  flight-path  angle  are  preserved  by 
transforming  the  state  vector  from  Earth  Centered  Inertial  to  the  Earth  Centered  Earth 
Fixed  reference  frame  as  defined  in  Appendix  C. 


Re-entry  Vehicle  Ground-level  Footprint  with  0°  Bank-angle 


Figure  11.  Range  of  a  re-entry  vehicle  with  and  without  the  rotation  of  the  Earth 

included  in  the  dynamics. 

4.3  Case  1  Results  and  Analysis 

The  optimal  solution  to  the  TPBVP  with  constraints  is  solved  using  the  GPM.  An 
initial  "point-to-target"  guess  is  provided  to  the  solver  by  integrating  the  initial 
conditions  forward  for  2000  seconds.  To  simplify  the  integration,  the  bank-angle  is  fixed 
to  zero  degrees.  This  choice  of  initial  guess  is  simply  an  engineering  judgment. 
Integrating  the  initial  conditions  is  a  means  of  providing  a  state  history  that  satisfies  the 
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equations  of  motion,  and  2000  seconds  (~33  min)  is  a  reasonable  estimate  of  the  flight 
duration. 

The  solution  obtained  from  the  first  iteration  of  the  GPM  solver  has  a  low 
number  of  nodes  (Recall,  it  is  easier  to  obtain  an  initial  solution  with  a  small  number  of 
nodes  compared  to  a  large  number  of  nodes),  and  as  discussed  earlier,  this  solution  is 
not  optimal.  The  solution  obtained,  however,  is  used  as  the  estimate  for  the  next 
iteration.  At  each  iteration  the  number  of  nodes  is  incrementally  increased  until  the 
minimum  of  sixty-one  nodes  is  reached.  Figure  12  shows  the  states  for  the  optimal 
solution  with  sixty-two  nodes.  In  the  last  subplot  showing  the  bank  angle,  the  RV 
appears  to  perform  a  rapid  bang-bang  maneuver  in  bank  angle  towards  the  end  of  the 
trajectory.  Investigation  into  this  behavior  shows  that  this  maneuver  is  much  smoother 
than  it  appears  (see  Section  4.6),  and  results  in  what  is  called  an  'S'  turn,  named  after  the 
apparent  path  of  the  vehicle  as  viewed  from  overhead.  This  last  minute  'S'  turn  allows 
the  RV  to  change  its  flight-path  angle  and  perform  a  dive  to  the  target  altitude.  The 

fifth  subplot  confirms  that  the  flight-path  angle  has  steepened,  reaching  nearly  -3° . 
Additionally,  the  dive  signifies  that  the  RV  is  dumping  energy  in  order  to  reach  the 
target  altitude,  but,  if  needed,  has  excess  energy  available  to  maneuver  to  additional 
waypoints.  The  use  of  this  excess  energy  to  reach  additional  waypoints  is  precisely  the 
motivation  behind  this  research. 
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Case  1:  Optimal  Solution  -  States  (62  nodes) 


800  1000  1200  1400  1600 


Figure  12.  Case  1  results  for  the  optimal  solution  without  predefined  waypoints. 


Once  the  baseline  solution  is  computed,  the  program  goes  about  the  task  of 
systematically  moving  the  waypoints  and  re-computing  a  solution.  The  GPM  solver 
automatically  interpolates  the  guess  to  fit  the  new  node  locations  using  cubic 
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interpolation.  Occasionally  the  interpolated  guess  results  in  exaggerated  states  (similar 
to  the  Runge  phenomenon),  and  the  GPM  solver  is  unable  to  converge  to  a  solution. 

One  method,  which  has  been  semi-successful,  is  to  select  another  set  of  nodes  to 
interpolate.  This  is  accomplished  in  the  code  by  incrementing  the  number  of  nodes  used 
in  the  solution.  Consequently  an  upper  node  limit  is  defined  so  that  the  algorithm  does 
not  continue  indefinitely  if  the  GPM  solver  is  continually  unsuccessful.  It  is  also 
expected  that  as  the  waypoints  are  pushed  further  away  from  the  nominal  trajectory,  an 
infeasible  solution  will  be  encountered.  An  infeasible  solution  is  one  in  which  the  RV 
does  not  have  enough  energy  or  maneuverability  to  both  reach  the  waypoint(s)  and  final 
target.  Thus,  the  upper  node  limit  also  serves  to  define  a  boundary  between  feasible  and 
infeasible  solutions. 

The  next  part  of  the  algorithm,  as  shown  in  Figure  7,  splits  the  optimal  solution 
into  two  phases  adjoined  at  the  latitude  and  longitude  corresponding  to  the  interior 
nodes  of  the  optimal  solution.  The  connection  of  the  two  phases  at  a  defined  latitude 
and  longitude  defines  a  waypoint.  Potentially  there  are  sixty-two  two-phase 
combinations;  however,  the  first  15%  and  last  15%  of  the  nodes  from  the  optimal 
solution  are  not  evaluated  under  the  assumption  that  the  variation  will  be  small  and  are 
of  less  interest  to  the  mission  planner.  Later  evidence  shows  that  this  assumption  is 
partially  true,  since  the  RV  is  still  capable  of  making  Targe'  maneuvers  near  the  end  of 
the  trajectory;  in  which  case,  the  last  15%  of  nodes  should  not  be  omitted.  For  this  study, 
of  the  original  sixty-two  nodes,  only  the  interior  forty-four  nodes  are  considered.  Each 
of  the  two-phase  combinations  is  then  re-solved  for  the  optimal  solution.  A  conclusion 
from  Bellman's  principle  of  optimality  translated  to  this  current  problem  states  that 
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starting  at  any  point  on  the  optimal  trajectory  results  in  the  new  optimal  trajectory 
following  the  same  optimal  trajectory  [8],  As  a  result,  there  are  ideally  forty-four 
identical  trajectories. 

The  algorithm  then  moves  the  waypoints  north  and  south  in  increments  of 
0.01°  (~0.6  nmi)  and  re-solves  the  optimal  control  problem.  The  results  from 
progressively  moving  the  waypoints  away  from  the  nominal  solution  are  presented  in 
Figure  13.  The  right  facing  triangle  is  the  RV  starting  position,  and  the  dashed  line  is  the 
optimal  trajectory  without  defined  waypoints.  The  upright  and  upside-down  triangles 
mark  both  the  longitudinal  meridian  that  the  waypoint  is  pushed  along  and  the 
sampling  bounds.  For  case  1,  an  iteration  limit  of  100  is  imposed,  and  therefore,  the  total 
angular  difference  between  the  two  triangles  is  2°  (~120  nmi).  Stated  earlier,  the 
algorithm  continues  to  push  away  from  the  nominal  trajectory  until  an  upper  node  limit 
is  reached  or  until  the  algorithm  reaches  an  iteration  limit.  This  termination  decision 
results  in  two  prominent  features  as  shown  in  Figure  13.  There  is  a  smooth  expansion 
region  followed  by  a  straight  leg  and  finally  a  'chunk'  of  data  missing.  The  smooth 
expansion  is  indicative  of  a  boundary  between  feasible  and  infeasible  solutions,  whereas 
the  straight  leg  section  is  where  the  algorithm  reached  the  iteration  limit.  The  figure  also 
shows  a  'chunk'  of  missing  data  that  is  thought  to  be  an  anomaly  created  by  the 
algorithm  terminating  too  soon.  It  is  also  noted  that  the  contour  information  provided  is 
limited  to  the  region  of  waypoint  sampling.  This  last  comment  will  become  more 
apparent  later  in  case  3. 
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Case  1:  Contour  of  Terminal  Time  [sec]  wrt  Lat/Lon 


Ion  [deg] 


Figure  13.  Case  1  contour  of  terminal  time  w.r.t.  Lat/Lon. 


Figure  13  shows  the  cost  of  the  dynamic  optimization  problem  for  each 
waypoint.  As  waypoints  are  chosen  off  the  baseline  optimal  trajectory,  the  cost  (time  to 
the  final  target)  increases.  Further,  the  terminal  time  is  more  sensitive  to  waypoint 
location  early  on  in  the  trajectory  as  opposed  to  later  in  the  trajectory.  This  figure  does 
not  adequately  capture  the  contour  behavior  towards  the  end  of  the  trajectory,  but  closer 
analysis  shows  the  cost  does  become  more  sensitive  to  waypoint  selection  near  the  end 
of  the  trajectory. 
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Case  1 :  Contour  of  Terminal  Time  [sec]  wrt  Lat/Lon 
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Figure  14.  Magnification  of  Figure  133  highlighting  'missing'  piece. 


Figure  14  identifies  two  locations  that  the  algorithm  prematurely  triggered  a 
'stop'  condition  as  defined  in  Section  3.4.  This  is  a  fact  verified  by  plotting  the 
trajectories  from  the  neighboring  slices  and  observing  their  passage  through  the  empty 
region.  This  result  is  a  byproduct  of  automating  the  process  and  does  not  account  for 
every  exit  scenario.  It  is  wholly  possible  to  go  back  into  the  code  and  tweak  the  step-size 
or  other  parameters  until  a  solution  is  found,  but  the  exact  reason  why  the  program 
terminated  prematurely  is  lost  in  the  automation  process.  Ensuring  a  solution  when 
using  an  automated  process  is  a  topic  for  future  research. 
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Figure  15  exhibits  the  same  features  as  Figure  13  except  that  the  velocity 
decreases  as  waypoints  are  chosen  away  from  the  unrestricted  optimal  solution. 


Case  1:  Contour  of  Terminal  Velocity  [km/s]  wrt  Lat/Lon 


Ion  [deg] 


Figure  15.  Case  1  contour  of  terminal  velocity  w.r.t.  Lat/Lon. 


In  between  each  vertical  triangle  pair  exists  optimal  trajectories  that  pass  through 
neighboring  latitudes  at  the  same  longitude.  Waypoints  that  share  the  same  longitude 
are  a  family  of  waypoints  that  make  up  a  slice  of  the  region  of  feasible  waypoints 
(meridian  slices).  By  plotting  the  family  of  waypoints  we  can  visualize  the  sensitivity  of 
the  optimal  trajectory  at  distinct  longitudinal  locations.  Figure  16  to  Figure  18  show  the 
terminal  states  of  three  such  families  of  waypoints.  In  each  figure,  only  the  terminal 
time  and  velocity  appear  to  carry  the  same  parabolic  shape  from  one  family  of 
waypoints  to  another.  The  physical  interpretation  and  potential  significance  of  the 
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discontinuous  breaks  in  the  flight-path  angle  and  bank  angle  are  not  yet  understood,  but 
the  results  are  provided. 


Terminal  States  w.r.t.  Waypoints  along  -35.2219°  Longitude 
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Figure  16.  Case  1  terminal  states  with  waypoints  along  longitude  -35.2219°£' 


At  approximately  26.4°  latitude,  the  RV  has  a  slight  bump  up  in  velocity  which 
appears  to  directly  correspond  to  the  steeper  flight-path  angle.  Logically  this  makes 
sense;  the  steeper  flight-path  angle  allows  the  RV  to  fall  faster  to  the  Earth,  and  therefore 
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it  picks  up  additional  speed.  Note  also,  the  heading  angle  appears  to  change  nearly 
linearly  with  waypoint  location. 


Terminal  States  w.r.t.  Waypoints  along  -13.2159°  Longitude 


latitude  [deg] 


Figure  17.  Case  1  terminal  states  with  waypoints  along  longitude  - 1 3.2 1 59  E 


The  increase  in  velocity  seen  towards  the  latter  portion  of  the  trajectory  (right 
end)  is  associated  with  the  RV  performing  a  last  minute  dive  to  the  target  altitude. 
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Terminal  States  w.r.t.  Waypoints  along  2.2809°  Longitude 


Figure  18.  Case  1  terminal  states  with  waypoints  along  longitude  2.2809°£' 


As  in  the  previous  figures,  but  now  more  pronounced  in  Figure  18,  the 
discontinuous  breaks  in  the  flight-path  angle  clearly  correspond  to  the  discontinuous 
breaks  in  the  velocity. 
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The  terminal  states  as  depicted  in  the  meridian  slices  (longitudinal  variations  in 
waypoint)  show  the  parabolic  nature  of  the  terminal  time  and  velocity  associated  with 
waypoints  off  the  baseline  trajectory.  Additionally,  the  change  in  heading  angle  is 
approximately  linear  with  waypoint  displacement.  Having  demonstrated  a  technique  to 
show  the  sensitivity  on  the  terminal  state  for  the  non-rotating  Earth,  this  technique  and 
results  are  now  expanded  upon  in  cases  2  and  3. 

4.4  Case  2  Results  and  Analysis 

Case  2  presents  a  vehicle  similar  to  case  1;  however,  case  2  differs  from  case  1  in: 
mass  of  the  vehicle,  inertial  velocity  and  flight-path  angle,  due  to  the  rotation  rate  of  the 
Earth  being  turned  'on.'  The  optimal  states  are  shown  in  Figure  19.  As  a  note  to 
furthering  the  discussion  on  rotation  effects,  we  see  that  the  two  trajectories,  case  1  and 
case  2,  differ  by  approximately  400  seconds  and  a  terminal  speed  of  2  km/  s.  Although  a 
smaller  portion  of  this  difference  is  due  to  the  difference  in  vehicle  mass,  the  primary 
contribution  is  due  to  the  difference  in  inertial  velocities. 

The  surrounding  area  0.5°  (~30  nmi)  north  and  south  of  the  optimal  nodes  are 
sampled  in  increments  of  0.01  (~0.6  nmi).  Forty-three  slices  or  families  of  waypoints 
composing  a  total  of  3,640  evaluations  of  the  surrounding  field  are  compiled  together  to 
give  the  contours  for  terminal  time  and  velocity.  Figure  22  to  Figure  24  are  each  a  family 
of  waypoints  depicting  their  respective  terminal  states.  As  in  case  1,  the  terminal  time 
and  velocity  show  strong  parabolic  shapes  throughout.  The  discussion  of  the  figures  for 
cases  2  and  3  are  very  similar  to  that  in  case  1.  Where  appropriate,  unique  features  will 
be  discussed,  but  otherwise  the  reader  is  referred  to  case  1  Results  and  Analysis. 
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Case  2:  Optimal  Solution  -  States  (61  nodes) 
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Figure  19.  Case  2  results  for  the  optimal  solution  without  predefined  waypoints. 


Similar  to  the  previous  case,  the  RV  performs  a  last  minute  bank  maneuver  to 
reduce  the  lift,  and  reach  the  target  altitude. 
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Case  2:  Contour  of  Terminal  Time  [sec]  wrt  Lat/Lon 
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Figure  20.  Case  2  contour  of  terminal  time  w.r.t.  Lat/Lon. 
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Case  2:  Contour  of  Terminal  Velocity  [km/s]  wrt  Lat/Lon 
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Figure  21.  Case  2  contour  of  terminal  velocity  w.r.t.  Lat/Lon. 
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Terminal  States  w.r.t.  Waypoints  along  -54.6387°  Longitude 


Figure  22.  Case  2  terminal  states  with  waypoints  along  longitude  -54.638 TE 
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Terminal  States  w.r.t.  Waypoints  along  -27.1483°  Longitude 


Figure  23.  Case  2  terminal  states  with  waypoints  along  longitude  -21  E 
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Terminal  States  w.r.t.  Waypoints  along  -5.2312°  Longitude 
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Figure  24.  Case  2  terminal  states  with  waypoints  along  longitude  —5.23 12° is 


Similar  to  case  1,  the  terminal  states  as  depicted  in  the  meridian  slices  show  the 
parabolic  nature  of  the  terminal  time  and  velocity  associated  with  waypoints  off  the 
baseline  trajectory.  In  contrast,  the  linearity  of  the  heading  angle  used  in  case  1  appears 
to  no  longer  exist.  A  mission  planner  may  use  Figure  24  to  gain  understanding  how 
changing  a  waypoint  affects  the  terminal  state.  In  particular,  the  time  and  velocity  are 
primary  observables,  whereas  the  other  states  are  likely  to  significantly  change  which 
selection  of  subsequent  waypoints. 
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4.5  Case  3  Results  and  Analysis 

Case  3  is  an  extension  of  case  2.  It  addresses  the  question  of  how  the  trajectory 
sensitivity  changes  as  waypoints  are  successively  added,  and  put  in  the  perspective  of 
the  mission  planner,  it  answers  the  question,  "Where  else  can  the  RV  go?"  In  particular, 
case  3  has  one  predefined  intermediate  waypoint  selected  from  the  data  collected  in  case 
2.  Figure  25  shows  the  selected  waypoint  slightly  north  of  the  nominal  solution. 
Conveniently,  the  waypoint  selected  is  a  waypoint  previously  solved  in  case  2;  therefore, 
the  state  history  for  the  optimal  solution  through  this  waypoint  is  already  obtained  and 
shown  in  Figure  26. 


Case  3:  Selection  of  Waypoint  from  Contour  ofTerminal  Time  [sec] 
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Figure  25.  Case  3  selection  of  fixed  waypoint. 


The  circle  shown  in  the  inset  plot  is  the  fixed  waypoint  (— 22. 07°  is  ,21.59°  N )  that 
all  the  subsequent  calculated  trajectories  will  pass  through. 
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Case  3:  Optimal  Solution  -  States  (50/50  Nodes) 
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Figure  26.  Case  3  results  for  the  optimal  solution  with  a  predefined  waypoint. 


The  optimal  trajectory  through  the  fixed  waypoint  is  calculated  with  fifty  nodes 
in  each  phase  of  this  two-phase  solution  represented  by  the  blue  and  red  circles. 
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The  algorithm  used  to  solve  cases  1  and  2  is  repeated  here  with  the  assumption 
that  the  waypoints  are  selected  consecutively  in  time,  so  only  the  area  following  the 
waypoint  is  considered  in  the  analysis.  The  results  are  presented  in  Figure  27  and 
Figure  28.  These  are  the  figures  that  would  be  most  useful  to  a  mission  planner  when 
the  set  of  waypoints  are  overlaid  thereby  showing  which  waypoints  are  immediately 
feasible  and  what  the  cost  is  associated  with  flying  through  a  particular  waypoint. 


Case  3:  Contour  of  Terminal  Time  [sec]  wrt  Lat/Lon 
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Figure  27.  Case  3  contour  of  terminal  time  w.r.t.  Lat/Lon. 
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Case  3:  Contour  of  Terminal  Velocity  [km/s]  wrt  Lat/Lon 
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Figure  28.  Case  3  contour  of  terminal  velocity  w.r.t.  Lat/Lon. 
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Figure  29.  Case  3  composite  of  all  calculated  trajectories. 


Figure  29  is  an  'aerial  view'  composite  of  all  the  trajectories.  It  depicts  a  bowtie 
looking  silhouette  with  the  trajectories  all  crossing  at  the  fixed  waypoint  (as  it  must). 
Note,  the  altitude  through  the  waypoint  is  not  defined,  so  a  'side  view'  of  the  trajectories 
through  the  intermediate  waypoint  would  show  the  RV  ranging  from  78  to  81  km  in 
altitude. 


63 


Terminal  States  w.r.t.  Waypoints  along  -12.3663°  Longitude 
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Figure  30.  Case  3  terminal  states  with  waypoints  along  longitude  —12.3663 °  E 


Like  in  the  previous  meridian  slices.  Figure  30  to  Figure  32  show  that  the  time 
and  velocity  plots  maintain  a  consistent  parabolic  shape.  However,  the  flight-path  angle 
and  bank  angle  plots  are  noticeably  dissimilar  from  each  other.  In  Figure  30,  the  flight- 
path  angle  changes  smoothly  but  has  a  peculiar  region  in  which  the  curve's  concavity 
changes.  In  Figure  31,  the  flight-path  angle  is  strictly  concave  up,  and  Figure  32  has 
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discontinuous  breaks.  Again,  the  interpretation  and  potential  significance  of  these 
curves  are  not  yet  understood. 


Terminal  States  w.r.t.  Waypoints  along  -0.91823°  Longitude 
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Figure  31.  Case  3  terminal  states  with  waypoints  along  longitude  -0.91 82°  E 
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Terminal  States  w.r.t.  Waypoints  along  7.3697°  Longitude 


Figure  32.  Case  3  terminal  states  with  waypoints  along  longitude  7.3697°  E 


The  figures  just  presented  are  extremely  useful  to  a  mission  planner.  The 
contour  plots.  Figure  27  and  Figure  28,  identify  regions  where  an  RV  can  go  and  show 
the  immediate  cost  associated  with  maneuvering  to  another  waypoint.  Figure  30  to 
Figure  32  focus  on  specific  meridians  that  provide  additional  visual  insight  into  trends 
in  the  data  that  are  otherwise  obscured  by  too  much  data.  In  the  selection  of  additional 
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waypoints,  the  mission  planner  should  consider  small  deviations  in  the  trajectory  that 
may  significantly  affect  the  terminal  state.  These  terminal  state  characteristics  are 
captured  in  the  meridian  slices  which  identify  highly  sensitive  states  through  the 
discontinuous  jumps  in  the  data,  e.g.  the  flight-path  angle  and  bank  angle  at  5.55° 
latitude  in  Figure  32. 

4.6  Finalizing  Trajectory 

This  section  makes  some  concluding  observations  applicable  to  any  computed 
trajectory.  For  this  analysis,  a  trajectory  was  selected  that  requires  the  RV  to  fly  through 
two  waypoints.  The  solution  obtained  is  presented  in  Figure  33  where  the  state  history 
is  shown  at  discrete  points  in  time.  At  approximately  1900  seconds,  it  appears  that  there 
is  a  discontinuous  break  in  the  bank  angle  data.  Flowever,  the  RV  could  rotate  as  fast  as 
a  =  ±60  °/sec ,  as  defined  in  Section  3.2.5.  The  spacing  of  the  nodes  around  1900 
seconds  prevents  an  adequate  capture  of  the  vehicle  dynamics.  A  solution  could  be  to 
simply  fit  a  cubic  polynomial  through  the  data,  propagate  the  initial  state  vector  using 
the  interpolated  control  data,  and  then  compare  to  see  if  the  two  solutions  match.  Or, 
the  problem  can  be  re-solved  using  more  nodes  so  that  the  vehicle  dynamics  are 
captured.  The  latter  option  has  the  benefit  of  improving  the  Hamiltonian,  and  therefore 
improving  the  overall  optimality  of  the  solution.  Figure  34  shows  the  states  again,  but 
this  time  using  four  phases  with  ninety  nodes  in  each.  At  approximately  1900  seconds, 
the  bank  angle  no  longer  appears  discontinuous.  We  also  see  in  Figure  35  that  the 
Hamiltonian  for  the  four  phase  solution  is  more  tightly  packed  to  the  optimal  H  =  —  1 . 
The  point  of  this  last  section  is  to  further  investigate  some  'features'  of  the  previous 
results.  A  numerical  solution  does  not  tell  the  whole  truth,  as  might  an  analytical 
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solution.  It  is  necessary  to  identify  at  what  point  in  the  trajectory  the  dynamics  are  not 
being  adequately  captured,  and  add  additional  nodes  to  fully  capture  the  dynamics.  A 
check  of  the  Hamiltonian  will  verify  that  an  optimal  solution  has  been  found.  Figure  34 
shows  some  additional  detail  of  the  dynamic  behavior  while  Figure  35  illustrates  the 
optimality  of  the  results.  This  section  builds  confidence  in  the  results  of  the  earlier 
sections. 


68 


Optimal  Trajectory  Through  2  Waypoints 
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Figure  33.  Optimal  trajectory  through  two  waypoints  represented  by  three  phases  (55- 

55-55  nodes). 
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Optimal  Trajectory  Through  2  Waypoints  -  360  nodes 
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Figure  34.  Optimal  trajectory  through  two  waypoints  using  four  phases  (90-90-90-90 

nodes) 


70 


Hamiltonian 


-0.95q 

I! 

-0  96Ci 
-0.97C3 
-0.98  - 


ti 


-0  994, 
-1 


-1  01  - 


-1.02  - 


O  3  phase  - 155  nodes 
□  4  phase  -  360  nodes 


1] 


o  cl 


0  200  400  600  800  1000  1200  1400  1600  1800  2000 

time  [sec] 


Figure  35.  Hamiltonian  of  solution  with  155  and  360  nodes. 
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4.7  Summary  of  Analysis 

An  incremental  approach  is  taken  to  determine  the  sensitivity  of  the  terminal 
state  to  waypoint  location.  A  minimum  number  of  nodes  used  in  the  solutions  is 
established  by  propagating  the  initial  conditions  through  an  integrator  and  increasing 
the  nodes  until  the  pseudospectral  program  output  coincides  with  the  integrated 
solution.  Comparing  the  footprints  of  straight,  constant  bank-angle  flight  trajectories 
with  the  Earth  rotation  'off'  and  'on'  justify  the  need  to  incorporate  the  Earth's  rotation 
rate  into  the  equations  of  motion.  After  establishing  an  optimal  trajectory  between  the 
initial  and  final  target,  similar  optimal  control  problems  are  solved  with  the  restriction 
that  the  RV  passes  through  a  nearby  waypoint.  Solving  the  optimal  control  problem  at 
numerous  neighboring  locations  and  tabulating  the  terminal  state  of  each  then  describes 
the  sensitivity  of  the  original  optimal  trajectory  to  changes  in  waypoint  location.  The 
process  was  fully  automated  to  select  the  neighboring  waypoints,  solve  the  optimal 
trajectory  for  each  waypoint,  and  then  construct  the  contour  plots.  The  changes  to  the 
final  state  as  a  function  of  waypoint  location  are  plotted  for  each  case. 
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V.  Conclusions  and  Future  Work 


5.1  Conclusions 

In  the  course  of  solving  the  optimal  control  problem,  this  research  built  upon  the 
already  successful  approach  using  GPM  presented  in  [23]  by  including  the  spherical 
Earth  model  and  the  rotation  rate  of  the  Earth.  With  the  addition  of  the  new  dynamics, 
it  was  realized  that  a  3D  flat,  non-rotating  Earth  is  not  suitable  for  accurately  modeling 
shallow  angle  orbital  re-entry  problems.  The  increased  model  fidelity,  however,  also 
means  abandoning  absolute  verification  of  the  trajectory's  optimality  as  analytical 
solutions  are  no  longer  available,  although  some  assurance  to  the  optimality  can  be 
gained  by  checking  the  Hamiltonian. 

The  introduction  of  uniformly  spaced  waypoints  'north'  and  'south'  of  the 
nominal  optimal  trajectory  serve  to  sample  the  surrounding  field.  Each  waypoint 
location  is  a  complete  optimal  control  problem  whereby  a  new  optimal  trajectory  is 
computed,  satisfying  the  waypoint  and  terminal  constraints.  Recognizing  that  the 
objective  cost  is  the  terminal  time,  comparing  the  terminal  states  of  the  baseline  solution 
to  the  waypoint  specified  trajectories  is  a  direct  measurement  of  the  system's  sensitivity 
to  waypoint  location.  The  process  was  fully  automated  so  that  the  user  does  not  need  to 
have  knowledge  of  optimal  control  or  need  to  manually  select  waypoints.  Even  with 
automation,  computation  time  for  problems  presented  in  this  research  (each  case) 
ranged  from  three  to  ten  days  on  a  Pentium  dual  core  32-bit  2  GHz  processor. 

5.2  Future  Work 

There  are  four  main  avenues  of  future  work:  continue  increasing  the  model 
fidelity  to  the  re-entry  vehicle  and  environment,  expand  the  sampling  region  to  include 
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all  feasible  solutions,  add  additional  features  for  specific  mission  requirements,  and  seek 
faster  solutions. 

Obvious  enhancements  to  model  fidelity  include  applying  the  heating  constraint, 
using  angle-of-attack  as  a  control,  and  CL  and  CD  as  a  function  of  angle-of-attack.  These 
three  modifications  would  round  out  the  basic  model  for  any  generic  re-entry  vehicle. 
Once  the  vehicle  is  well-defined,  other  considerations  such  as  using  better  atmospheric 
and  Earth  models  are  within  reason.  Beyond  these  suggestions  is  an  exercise  in 
applying  corrections  on  top  of  corrections. 

It  is  assumed  that  the  range  of  all  feasible  solutions  is  continuous  within  an 
enclosed  envelope  surrounding  the  nominal  trajectory.  This  research,  although  not 
focused  on  finding  the  extreme  upper  and  lower  boundaries,  has  provisions  to  continue 
until  unable  to  continue  further.  At  each  slice,  the  algorithm  is  instructed  to  stop 
proceeding  after  a  maximum  number  of  increments  or  an  infeasible  solution  is  reached. 
If  the  maximum  iteration  limit  restriction  is  lifted,  an  upper  and  lower  bound  of  feasible 
solutions  along  each  slice  is,  in  theory,  obtainable.  The  'stop'  conditions  mentioned 
earlier  resulted  in  the  program  terminating  too  soon.  Aside  from  the  normal  'stop' 
conditions,  the  event  may  be  triggered  by  too  large  a  step  size,  a  poor  initial  guess  at  the 
solution,  or  a  hypothetical  infeasible  region  contained  within  the  envelope  of  all  feasible 
solutions.  Results  presented  in  Chapter  IV  suggest  such  infeasible  regions  may  exist; 
therefore,  future  work  in  identifying  the  region  of  all  feasible  waypoints  should  also 
answer  the  question  about  infeasible  regions. 

A  mission  feature  omitted  from  this  research  but  included  in  Maj  Jorris' 
dissertation  is  the  presence  of  no-fly  zones  built  directly  into  the  optimization  problem. 
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Therefore  the  reincorporation  of  no-fly  zones  and  previously  recommended  future  work 
from  [23]  are  certainly  viable. 

The  initial  algorithm  and  design  for  this  research  was  laid  out  early  on  under  the 
assumption  that  sampling  the  field  at  many  points  would  not  take  very  long  to  compute. 
This  decision  was  due  to  expectations  set  by  the  method  converging  to  an  optimal 
solution  in  seconds  from  the  work  conducted  in  [23],  The  algorithm  presented  in  Figure 
7  is  constructed  in  such  a  way  that  adaptation  to  a  distributed  or  parallel  computing 
system  is  feasible.  Obtaining  results  quicker  than  the  last  is  always  heralded  as  a  good 
achievement,  and  so  future  research  should  seek  to  reduce  the  computation  time  (in  lieu 
of  processor  speed  advancements)  by  means  such  as  algorithmic  changes,  problem 
reformulation,  and  efficient  coding. 
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Appendix  A.  Non-dimensionalization 


A.1  Non-dimensional  3-D  Re-entry  Equations  of  Motion  -  Rotating  Earth 


The  3D  equations  of  motion  from  [20],  omitting  the  thrust  terms,  are: 
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Note  gravity  is  a  function  of  r,  and  can  be  rewritten  in  terms  of  the  new  non-dimensional 
term  r  . 
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Also  the  rotation  rate  of  the  Earth,  C0m ,  needs  to  be  non-dimensionalized  which  can  be 


achieved  by  multiplying  by  J  y  .  Therefore  we  now  have: 
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Now  the  non-dimensional  equations  of  motion  can  be  derived: 
Radial  Distance: 
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Velocity: 
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1 

^3 

(A 

i'o/ 

/ Sq  j 

■\fa/o  V 

O 

1 

0 

^0 

(A.  17) 


dRV'  _  1 
dr  g0 


- g  sin  y  +  rO)l  cos  <p  ( cos  (/>  sin  y  -  sin  <j)  sin  y  cos  y ) 

m 


(A.18) 


dRV' 

dr 


sin  y  |  :  /  :  .  ,  .  \  J--./ 

— 77-  +  r  cos  <p  ( cos  (p  sin  y  -  sm  (p  sin  y  cos  y)-D 

r " 


(A.19) 


Longitude: 


Latitude: 


dG'  dO 


/  ro  dd 

yg  0  dt 


dO' 

dr 


RV  cos  y  cos  y 
r  cos  (p 


da'  _ 

r  Ry'^ 

cosy  cos  y 

dr 

l  r) 

cos^ 

d£_ 

dr 


d£ 

dr 


d<p' 

r  Ryy 

-  — 

cosy  sin  y 

dr 

l  r  J 

(A.20) 

(A.21) 

(A.22) 

(A.23) 

(A.24) 

(A.25) 
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Flight-path-angle: 


dy 

dr 


dy 

dr 


g  V  _  , 

—cos  y-\ - cos^+2<y0  cos^cos^ 

V  r 


+  ... 


...+ 


raj 

RV 


cos  (j)[  cos  </>  cos  sin  (j)  sin  y/  sin  y) 


(A.26) 


(A.27) 


dr  L'  cos  y  RV'  .  , 

—  =  — COS<7 — 7^-7  +  —  COS^+2tye  cos^cos^  +  ... 
dT  V  r  V  r 

/ 

V 

. . .  +  &42  cos  (j)  (cos  (j)  cos  y  +  sin  (/>  sin  y/  sin  y) 


(A.28) 


Heading: 


dy/' 

dr 


dy/' 

dr 


Lsinrr 
RVm  cos  y 


—  cos  ^cos  y/  tan  (f)  +  2co®  (sin  y/ cos ^  tan  y-  sin (/)) 
r 


r*b~  r°ji 

V  8o  Ry  cos  y 


sin  (j)  cos  (j)  cos  y/ 


(A.29) 


(A.30) 


dy/'  L'  sin  <t  RV'  ,  .  ,  ,  .  ,  .  , 

- =  - —  cos  ^cos  y/  tan  (p  +  2 co®  ( sin  y/  cos  (p  tan  y-  sin  (p ) . . . 

dr  V  cos  y  r 


r'cd®  sin  (j)  cos  (j)  cos  y/ 


(A.31) 


RV' 


cos  y 


The  atmosphere  is  modeled  as  a  simple  exponential. 


P  =  Pi)e 


~P(r~r0) 


(A.32) 


And  therefore  lift  can  be  rewritten 
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(A.33) 


L  = 


_  pCLSRV 2  _  p0e~^r~'°'1  C LS  RV2 


L  = 


_  pQe~^r'r°-r°'lCLS RV'2g0r0  _  p0e~Mr'~l)CLSRV'2gQr() 


And  dividing  by  g0m 


Similarly  drag  is: 


L'  = 


_  rQ 


Poe 


-Ao('-'-i) 


CLSRV'2gQr0 


8  om 


2  g0m 


L'  = 


Pi/' 


-Ao('-'-i) 


ClSrV'\ 


2m 


2/77 


(A.34) 


(A.35) 


(A.36) 


(A.37) 


In  summary  the  non-dimensional  equations  of  motion  are  (A.  16), (A.  19), (A. 22), (A. 25), 
(A.28),(A.31)  respectively: 
dr' 


R-it' 


V  sin  r 
dr 

dRV'  sin  y  ,2  ,,  ,  .  ...  A  _, 

- = - —  +  r  <y0" cos^(cos^sin ^-sin^sin^cos y)  —  D 

dr  r 


dO' 

dr 

d£ 

dr 


(  Ry'\ 


cosy  cos  y/ 


V  r  J 

f  Ry>\ 


V  r  J 


cos^ 

cos  y  sin  y/ 


ft  y-  /  Rt  j'  / 

dy  L  cos  y  V  „  ,  ,  r  ,2  ,,  .  ,  .  , 

- =  -^—7  cos  rr  — ,,,  -  ,  -1 — —  cos  y  +  2<w0  cos  <p  cos  iff  +  -^—7  ty0  cos  (p  ( cos  (p  cos  7  +  sin  <p  sin  yt  sin  y) 


dr  V  r  V  r  V 


dy/'  L'sinc  RV'  ,  _  ,  ,  .  ,  .  ,N  /<wl2  sin^cos^cost^ 

— —  =  -^—7 - T-cos7cos^tan^  +  2tyffi  (sin^cos^tan^-sin^)-- 

dr  V  cos  7  r 


cos  7 
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Appendix  B.  MATLAB®  Function  Descriptions 

Several  MATLAB®  scripts  and  functions  are  used  in  the  course  of  this  research. 
The  files  are  categorized  into:  problem  setup,  files  required  by  OPENPOCS,  and 
graphics  output.  A  MATLAB®  function  has  the  form: 

[  output  ]  =  function_name  ( input ) 

and  can  accept  multiple  inputs  and  outputs.  A  script  file  unlike  a  function  does  not 
require  any  input,  but  all  of  the  executed  code  contained  in  the  script  is  available  to  the 
calling  function. 

B.l  Problem  Setup 


main.m 

-  design_problem.m 

-  combine.m 

-  findsoln_ic2fc.m 

L  workhorse.m 

-  design_parameters.m 
L  <OPENPOCS> 

Hmdsoln_ic2wp2fc.m 

-  workhorse2.m 

-design_parameters.m 
L  <OPENPOCS> 
*-workhorse3.m 

-  design_parameters.m 
L  <OPENPOCS> 


Figure  36.  Program  directory  tree. 


•  main.m:  This  is  a  script  file  that  initializes  and  ties  together  all  the  subfunctions 
used  in  calculating  the  optimal  trajectories. 

•  design_problem.m:  This  script  file  contains  the  initial,  intermediate  (optional), 
and  final  conditions. 
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•  combine.m:  This  function  file  receives  as  input  the  initial,  intermediate 
(optional),  and  final  conditions  defined  in  desig7t_problem.nl,  and  combines  them 
into  a  single  structure. 

•  findsoln_ic2fc.m:  This  function  computes  the  unrestricted  optimal  solution 
between  the  initial  and  final  target. 

•  workhorse.m:  This  function  calculates  the  one-phase  optimal  trajectory  between 
the  initial  and  final  target. 

•  design_parameters.m:  This  is  a  script  file  that  contains  constants  for  the  model, 
vehicle,  and  scaling  factors. 

•  findsoln_ic2wp2fc.m:  This  function  splits  the  one-phase  optimal  solution 
obtained  in  findsoln_ic2fc.m  into  two-phase  solutions  and  selects  waypoints 
used  in  mapping  the  sensitivity  contours. 

•  workhorse2.m:  This  function  calculates  the  two-phase  optimal  trajectory 
between  the  initial  and  final  target  with  an  intermediate  waypoint  specified. 

•  workhorse3.m:  This  function  calculates  the  three-phase  optimal  trajectory 
between  the  initial  and  final  target  with  two  intermediate  waypoints  specified. 

B.2  Files  Required  By  OPENPOCS 

•  cost_mintime_openpocs.m:  This  function  contains  the  objective  cost  function. 

•  eom_hicks_nd_openpocs.m:  This  function  contains  the  equations  of  motion 
and  path  constraints. 

•  connection_openpocs.m:  This  function  specifies  how  the  phases  link  to 
eachother. 
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B.3  Graphics  Output 

•  printsoln.m:  This  function  prints  the  state  history  of  a  specific  trajectory. 

•  printcontours.m:  This  function  prints  the  terminal  time  and  velocity  contours. 
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Appendix  C.  Coordinate  Transformations 

C.l  Earth  Centered  Inertial  to  Earth  Centered  Earth  Fixed 

The  following  are  the  steps  required  to  convert  from  the  Earth  Centered  Inertial 
(ECI)  to  the  Earth  Centered  Earth  Fixed  Frame  (ECEF)  using  the  following  inertial  state 
vector: 


V 

e 

</> 

7 

Y_ 

Convert  position  in  inertial  frame  from  spherical  coordinates  to  the  Cartesian 
coordinates. 


(C.l) 


X 

rcos(y)cos(#) 

y 

= 

rcos(y)sin(#) 

1 

C/5 

_ 1 

(C.2) 


Likewise,  convert  the  velocity  from  spherical  to  Cartesian,  i.e.  the  east-north-up  frame. 


rv.i 

V  cos  (/)  cos  (y/) 

K 

= 

V  cos(^)  sin(^) 

y„\ 

V  sin  ( y) 

(C.3) 


Then  convert  the  velocity  to  the  inertial  frame.  This  is  accomplished  by  a  Body  Three  2-3 
rotation.  (Note:  the  coordinates  e-n-u  are  re-arranged  to  u-e-n) 


cos(-#) 

sin(-6>)  0 

cos(^)  0  -sin(^) 

K 

Vy 

= 

-sin(-0)  cos(-£?)  0 

0  1  0 

K 

yz_ 

0 

0  1 

sin(^)  0  cos(^) 

K_ 

The  velocity  due  to  Earth's  rotation  is  subtracted. 
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XI 

XI 

0 

0 

X 

v 

= 

X 

- 

0 

0 

0 

y 

(C.5) 

xJ 

XJ 

— 

0 

0 

z 

Convert  the  velocity  back  to  the  up-east-north  frame. 


XI 

cos(-^)  0  -sin(-0) 

cos(#) 

sin(#)  0 

X" 

v. 

= 

0  1  0 

-sin(#)  cos(#)  0 

X 

XJ 

sin(-^)  0  cos(-^) 

0 

0  1 

X. 

(C.6) 


And  finally  convert  the  velocity  back  to  spherical  coordinates,  taking  care  to  recognize  a 
quadrant  check  is  required. 


>1 

tan -'(VJV„) 

X 

= 

tan -'(vJfiT+vJ 
1  - 

V 

A. 

The  relative  state  vector  is  now: 


(C.7) 


V 

e 

Rr 

¥ 


(C.8) 
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